IRAP trial-type D scores are calculated from an average of only 18 pairs of reaction times. This would be deemed as far too low anywhere else in the literature on reaction time based tasks. The implications of this can be seen in how poorly estimated any one IRAP D score is. We can observe this by bootstrapping reaction times for each participant’s D scores and noting how wide their confidence intervals are.
“Credibility Interval” is used in a number of different ways in the literature, sometimes denoting the difference between frequentist vs Bayesian approaches to probability, and other times denoting whether the interval takes into account not only the variance of the fixed effect in question but also the random effects. Here I employ the latter sense of the word.
The metafor R package easily produces credibility intervals but typically employed effect sizes as its input, whereas here I employ lme4 to fit mutilevel/metaanalytic models using the data itself. lme4 and helper packages built on top of it do not provide a method to produce credibility intervals in the sense that metafor does. As such, I implemented this myself using the metafor package’s code. The key conceptual link between the two is that lme4 models return the SD of the random effect, whereas metafor refers to this as Tau and return it as its squared value. These values refer to the same property: “In random-effects meta-analysis, the extent of variation among the effects observed in different studies (between-study variance) is referred to as tau-squared, τ2, or Tau2 (Deeks et al 2008). τ2 is the variance of the effect size parameters across the population of studies and it reflects the variance of the true effect sizes. The square root of this number is referred to as tau (T). T2 and Tau reflect the amount of true heterogeneity. T2 represents the absolute value of the true variance (heterogeneity). T2 is the variance of the true effects while tau (T) is the estimated standard deviation of underlying true effects across studies (Deeks et al 2008). The summary meta-analysis effect and T as standard deviation may be reported in random-effects meta-analysis to describe the distribution of true effects (Borenstein et al 2009).” (Source).
# dependencies
library(tidyverse)
library(knitr)
library(kableExtra)
library(boot)
library(parallel)
library(bayestestR)
library(patchwork)
library(mdthemes)
library(lme4)
library(sjPlot)
library(emmeans)
library(ggstance)
# library(merTools) called via merTools:: to avoid namespace collisions between MASS and dplyr
# set seed for reproducibility
set.seed(42)
# number of bootstraps
n_boots <- 2000
# options
options(knitr.table.format = "html") # necessary configuration of tables
# disable scientific notation
options(scipen = 999)
# function to round all numeric vars in a data frame
round_df <- function(df, n_digits = 3) {
df %>% mutate_if(is.numeric, round, digits = n_digits)
}
# get data from evaluative IRAPs
data_trial_level <- read_csv("../data/data_trial_level.csv") %>%
filter(timepoint == "baseline")
# outliers
data_outliers <- data_trial_level %>%
distinct(unique_id, .keep_all = TRUE) %>%
select(unique_id, domain, mean_rt) %>%
mutate(median_mean_rt = median(mean_rt, na.rm = TRUE),
mad_mean_rt = mad(mean_rt, na.rm = TRUE)) %>%
# exclude median +- 2MAD
mutate(rt_outlier = ifelse(mean_rt < median_mean_rt-mad_mean_rt*2 |
mean_rt > median_mean_rt+mad_mean_rt*2, TRUE, FALSE)) %>%
filter(rt_outlier == FALSE) %>%
select(unique_id, rt_outlier) %>%
full_join(data_trial_level, by = "unique_id") %>%
mutate(rt_outlier = ifelse(is.na(rt_outlier), TRUE, rt_outlier))
data_outliers_removed <- data_outliers %>%
filter(rt_outlier == FALSE)
# trim RTs>10000 ms, as part of D scoring
data_trimmed <- data_outliers_removed %>%
select(unique_id, domain, trial_type, rt, block_type) %>%
filter(rt <= 10000)
# for plot
dodge_width <- 0.25data_outliers %>%
distinct(unique_id, .keep_all = TRUE) %>%
count(rt_outlier) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| rt_outlier | n |
|---|---|
| FALSE | 1464 |
| TRUE | 110 |
data_descriptives <- data_outliers_removed %>%
distinct(unique_id, .keep_all = TRUE)
data_descriptives %>%
count(domain) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| domain | n |
|---|---|
| Body image | 23 |
| Clinton-Trump | 97 |
| Countries (1) | 53 |
| Countries (2) | 50 |
| Death (1) | 17 |
| Death (2) | 20 |
| Death (3) | 26 |
| Disgust (1) | 37 |
| Disgust (2) | 43 |
| Friend-Enemy | 98 |
| Gender (1) | 41 |
| Gender (2) | 95 |
| Lincoln-Hitler | 131 |
| Personality - Agreeableness | 39 |
| Personality - Conscientiousness | 39 |
| Personality - Extraversion | 33 |
| Personality - Neuroticism | 32 |
| Personality - Openness | 33 |
| Race (1) | 45 |
| Race (2) | 44 |
| Religion | 30 |
| Rich-Poor | 84 |
| Sexuality (1) | 26 |
| Sexuality (2) | 19 |
| Shapes & colors (1) | 11 |
| Shapes & colors (2) | 11 |
| Shapes & colors (3) | 10 |
| Shapes & colors (4) | 8 |
| Shapes & colors (5) | 22 |
| Shapes & colors (6) | 26 |
| Shapes & colors (7) | 14 |
| Stigma - ADHD | 62 |
| Stigma - PTSD | 54 |
| Stigma - Schizophrenia | 54 |
| Valenced words | 37 |
data_descriptives %>%
count(domain) %>%
summarize(total_n = sum(n),
min_n_per_domain = min(n),
max_n_per_domain = max(n),
mean_n_per_domain = round(mean(n, na.rm = TRUE), 2),
sd_n_per_domain = round(sd(n, na.rm = TRUE), 2)) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| total_n | min_n_per_domain | max_n_per_domain | mean_n_per_domain | sd_n_per_domain |
|---|---|---|---|---|
| 1464 | 8 | 131 | 41.83 | 28.73 |
data_descriptives %>%
summarize(min_age = round(min(age, na.rm = TRUE), 2),
max_age = round(max(age, na.rm = TRUE), 2),
mean_age = round(mean(age, na.rm = TRUE), 2),
sd_age = round(sd(age, na.rm = TRUE), 2)) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| min_age | max_age | mean_age | sd_age |
|---|---|---|---|
| 16 | 60 | 20.1 | 4.39 |
data_descriptives %>%
count(gender) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| gender | n |
|---|---|
| Female | 221 |
| Male | 134 |
| Other | 1 |
| NA | 1108 |
ggplot(data_trimmed, aes(rt, fill = block_type)) +
geom_density(alpha = 0.3) +
facet_wrap(~ trial_type)data_trimmed %>%
group_by(trial_type, block_type) %>%
summarize(mean = mean(rt, na.rm = TRUE),
sd = sd(rt, na.rm = TRUE)) %>%
round_df(0) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| trial_type | block_type | mean | sd |
|---|---|---|---|
| tt1 | con | 1395 | 567 |
| tt1 | incon | 1560 | 637 |
| tt2 | con | 1570 | 658 |
| tt2 | incon | 1633 | 687 |
| tt3 | con | 1560 | 639 |
| tt3 | incon | 1520 | 670 |
| tt4 | con | 1598 | 662 |
| tt4 | incon | 1634 | 695 |
# bootstrapping has a long execution time, so load saved values if they've already been calculated
if(file.exists("models/data_estimates_D.rds")) {
data_estimates_D <- read_rds("models/data_estimates_D.rds")
} else {
D_score <- function(data, i) {
data_with_indexes <- data[i,] # boot function requires data and index
mean_con <- mean(data_with_indexes$rt[data_with_indexes$block_type == "con"], na.rm = TRUE)
mean_incon <- mean(data_with_indexes$rt[data_with_indexes$block_type == "incon"], na.rm = TRUE)
sd <- sd(data_with_indexes$rt, na.rm = TRUE)
D <- (mean_incon - mean_con) / sd
return(D)
}
bootstrap_D_score <- function(data){
require(dplyr)
require(boot)
fit <-
boot::boot(data = data,
statistic = D_score,
R = n_boots,
sim = "ordinary",
stype = "i",
parallel = "multicore",
ncpus = parallel::detectCores())
results <- boot::boot.ci(fit, conf = 0.95, type = c("norm","basic", "perc", "bca"))
output <-
tibble(method = c("normal", "basic", "percent", "bca"),
estimate = rep(fit$t0, 4),
ci_lower = c(results$normal[2], results$basic[4], results$percent[4], results$bca[4]),
ci_upper = c(results$normal[3], results$basic[5], results$percent[5], results$bca[5]))
return(output)
}
# bootstrap D scores
data_estimates_D <- data_trimmed %>%
select(unique_id, domain, trial_type, rt, block_type) %>%
group_by(unique_id, domain, trial_type) %>%
do(bootstrap_D_score(data = .)) %>%
ungroup() %>%
mutate(sig = ifelse((ci_lower < 0 & ci_upper < 0) | (ci_lower > 0 & ci_upper > 0), TRUE, FALSE),
ci_width = ci_upper - ci_lower) %>%
round_df(3)
# save to disk
write_rds(data_estimates_D, "models/data_estimates_D.rds", compress = "gz")
}# bootstrapping has a long execution time, so load saved values if they've already been calculated
if(file.exists("models/data_estimates_PI.rds")) {
data_estimates_PI <- read_rds("models/data_estimates_PI.rds")
} else {
# Fast calculation of the A statistic - code from Ruscio (2008) supplementary materials
PI_score <- function(data, i) {
data_with_indexes <- data[i,] # boot function requires data and index
x <- na.omit(data_with_indexes$rt[data_with_indexes$block_type == "incon"])
y <- na.omit(data_with_indexes$rt[data_with_indexes$block_type == "con"])
nx <- length(x)
ny <- length(y)
rx <- sum(rank(c(x, y))[1:nx])
PI <- (rx / nx - (nx + 1) / 2) / ny
return(PI)
}
bootstrap_PI_score <- function(data){
require(dplyr)
require(boot)
fit <-
boot::boot(data = data,
statistic = PI_score,
R = n_boots,
sim = "ordinary",
stype = "i",
parallel = "multicore",
ncpus = parallel::detectCores())
results <- boot::boot.ci(fit, conf = 0.95, type = c("norm","basic", "perc", "bca"))
output <-
tibble(method = c("normal", "basic", "percent", "bca"),
estimate = rep(fit$t0, 4),
ci_lower = c(results$normal[2], results$basic[4], results$percent[4], results$bca[4]),
ci_upper = c(results$normal[3], results$basic[5], results$percent[5], results$bca[5]))
return(output)
}
# bootstrap PI scores
data_estimates_PI <- data_outliers_removed %>%
group_by(unique_id, domain, trial_type) %>%
do(bootstrap_PI_score(data = .)) %>%
ungroup() %>%
mutate(sig = ifelse((ci_lower < 0.50 & ci_upper < 0.50) | (ci_lower > 0.50 & ci_upper > 0.50), TRUE, FALSE),
ci_width = ci_upper - ci_lower) %>%
round_df(3)
# save to disk
write_rds(data_estimates_PI, "models/data_estimates_PI.rds", compress = "gz")
}data_estimates_D %>%
arrange(estimate) %>%
mutate(method = fct_relevel(method, "bca", "basic", "normal", "percentile")) %>%
group_by(method) %>%
mutate(ordered_id = row_number()/n()) %>%
ungroup() %>%
ggplot() +
geom_linerange(aes(x = ordered_id, ymin = ci_lower, ymax = ci_upper, color = sig),
alpha = 1) +
geom_point(aes(ordered_id, estimate), size = 0.5) +
geom_hline(yintercept = 0, linetype = "dotted") +
mdthemes::md_theme_linedraw() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "top") +
scale_color_viridis_d(end = 0.6, direction = -1) +
xlab("Participant (ranked by *D* score)") +
ylab("*D* score") +
labs(color = "95% CI excludes zero point") +
facet_wrap(~ method, ncol = 5)p_cis_by_domain <-
data_estimates_D %>%
filter(method == "bca") %>%
mutate(domain = str_replace(domain, "Personality - ", "Big5: "),
domain = str_replace(domain, "Stigma - ", "Stigma: ")) %>%
arrange(estimate) %>%
group_by(domain) %>%
mutate(ordered_id = row_number()/n()) %>%
ungroup() %>%
ggplot() +
geom_linerange(aes(x = ordered_id, ymin = ci_lower, ymax = ci_upper, color = sig),
alpha = 1) +
geom_point(aes(ordered_id, estimate), size = 0.5, shape = "square") +
geom_hline(yintercept = 0, linetype = "dotted") +
mdthemes::md_theme_linedraw() +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "top") +
scale_color_viridis_d(end = 0.6, direction = -1) +
xlab("Ranked participant") +
ylab("*D* score") +
labs(color = "95% CI excludes zero point") +
facet_wrap(~ domain, ncol = 5)
p_cis_by_domainWidths cant be directly compared between D and PI as they have different ranges, so D scores only.
Not meta analyzed as extreme skew in data means that residuals are very non-normal, violating assumptions and underestimating MAP estimates. Instead I simply present MAP estimates for each trial type & method, and then domain, trial type, and method.
full_join(data_estimates_D %>%
mutate(method = fct_relevel(method, "bca", "basic", "normal", "percent")) %>%
group_by(method) %>%
do(point_estimate(.$ci_width, centrality = "MAP")),
data_estimates_D %>%
mutate(method = fct_relevel(method, "bca", "basic", "normal", "percent")) %>%
group_by(method) %>%
summarize(min = min(ci_width),
max = max(ci_width),
.groups = "drop"),
by = "method") %>%
round_df(2) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| method | MAP | min | max |
|---|---|---|---|
| bca | 1.31 | 0.34 | 1.97 |
| basic | 1.31 | 0.39 | 1.75 |
| normal | 1.32 | 0.37 | 2.15 |
| percent | 1.31 | 0.39 | 1.75 |
data_estimates_D %>%
mutate(method = fct_relevel(method, "bca", "basic", "normal", "percent")) %>%
group_by(trial_type, method) %>%
do(point_estimate(.$ci_width, centrality = "MAP")) %>%
round_df(2) %>%
pivot_wider(names_from = method, values_from = MAP) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| trial_type | bca | basic | normal | percent |
|---|---|---|---|---|
| tt1 | 1.31 | 1.31 | 1.32 | 1.31 |
| tt2 | 1.31 | 1.31 | 1.32 | 1.31 |
| tt3 | 1.31 | 1.31 | 1.32 | 1.31 |
| tt4 | 1.31 | 1.31 | 1.32 | 1.31 |
data_estimates_D %>%
mutate(method = fct_relevel(method, "bca", "basic", "normal", "percent")) %>%
summarize(min_Dscore = min(estimate),
max_Dscore = max(estimate)) %>%
mutate(range_Dscore = max_Dscore - min_Dscore) %>%
round_df(2) ## # A tibble: 1 x 3
## min_Dscore max_Dscore range_Dscore
## <dbl> <dbl> <dbl>
## 1 -1.27 1.53 2.8
data_ci_width_map_D <- data_estimates_D %>%
group_by(domain, trial_type, method) %>%
#summarize(median = median(ci_width), .groups = "drop") %>%
do(point_estimate(.$ci_width, centrality = "MAP")) %>%
ungroup() %>%
mutate(method = fct_relevel(method, "bca", "basic", "normal", "percent"),
method = fct_rev(method),
trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
trial_type == "tt2" ~ "Trial type 2",
trial_type == "tt3" ~ "Trial type 3",
trial_type == "tt4" ~ "Trial type 4"),
trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4")) %>%
mutate(domain = fct_reorder(domain, MAP, .fun = min)) %>%
arrange(domain)
# save to disk
write_rds(data_ci_width_map_D, "models/data_ci_width_map_D.rds")
# plot
p_ci_widths <-
ggplot(data_ci_width_map_D, aes(MAP, domain, color = method, shape = method)) +
geom_line() +
geom_point() +
scale_shape_manual(labels = c("BCA", "Basic", "Normal", "Percentile"), values = c(15, 16, 17, 7)) +
scale_color_viridis_d(begin = 0.2, end = 0.8, labels = c("BCA", "Basic", "Normal", "Percentile")) +
mdthemes::md_theme_linedraw() +
facet_wrap(~ trial_type, ncol = 4, nrow = 1) +
labs(x = "MAP 95% CI width",
y = "",
color = "Bootstrap method",
shape = "Bootstrap method") +
theme(legend.position = "top")
p_ci_widthsdata_diff_zero <-
bind_rows(mutate(data_estimates_D, DV_type = "*D* scores"),
mutate(data_estimates_PI, DV_type = "PI scores")) %>%
mutate(domain = as.factor(domain),
method = fct_relevel(method, "bca", "basic", "normal", "percent"),
method = fct_rev(method),
trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
trial_type == "tt2" ~ "Trial type 2",
trial_type == "tt3" ~ "Trial type 3",
trial_type == "tt4" ~ "Trial type 4"),
trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4")) %>%
group_by(domain, trial_type, method, DV_type) %>%
summarize(proportion_diff_zero = mean(sig),
variance = plotrix::std.error(sig)^2,
.groups = "drop") %>%
mutate(variance = ifelse(variance == 0, 0.000001, variance)) # model cannot be run on zero variance, so offset by a miniscule amount
# save to disk
write_rds(data_diff_zero, "models/data_diff_zero.rds")data_diff_zero %>%
filter(method == "bca") %>%
mutate(DV_type = case_when(DV_type == "*D* scores" ~ "D",
DV_type == "PI scores" ~ "PI")) %>%
select(-variance) %>%
pivot_wider(names_from = c(trial_type, DV_type, method),
values_from = proportion_diff_zero) %>%
round_df(2) %>%
janitor::clean_names() %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| domain | trial_type_1_d_bca | trial_type_1_pi_bca | trial_type_2_d_bca | trial_type_2_pi_bca | trial_type_3_d_bca | trial_type_3_pi_bca | trial_type_4_d_bca | trial_type_4_pi_bca |
|---|---|---|---|---|---|---|---|---|
| Body image | 0.17 | 0.17 | 0.09 | 0.09 | 0.04 | 0.13 | 0.22 | 0.13 |
| Clinton-Trump | 0.20 | 0.25 | 0.08 | 0.07 | 0.18 | 0.18 | 0.08 | 0.05 |
| Countries (1) | 0.21 | 0.34 | 0.15 | 0.19 | 0.13 | 0.08 | 0.00 | 0.02 |
| Countries (2) | 0.08 | 0.14 | 0.12 | 0.20 | 0.08 | 0.12 | 0.06 | 0.06 |
| Death (1) | 0.06 | 0.06 | 0.12 | 0.12 | 0.29 | 0.18 | 0.06 | 0.06 |
| Death (2) | 0.10 | 0.10 | 0.10 | 0.10 | 0.10 | 0.15 | 0.10 | 0.20 |
| Death (3) | 0.38 | 0.31 | 0.15 | 0.19 | 0.12 | 0.08 | 0.15 | 0.19 |
| Disgust (1) | 0.19 | 0.24 | 0.11 | 0.11 | 0.30 | 0.24 | 0.32 | 0.30 |
| Disgust (2) | 0.23 | 0.26 | 0.21 | 0.23 | 0.12 | 0.07 | 0.33 | 0.35 |
| Friend-Enemy | 0.31 | 0.36 | 0.18 | 0.18 | 0.04 | 0.09 | 0.08 | 0.11 |
| Gender (1) | 0.07 | 0.10 | 0.00 | 0.05 | 0.10 | 0.07 | 0.17 | 0.20 |
| Gender (2) | 0.19 | 0.17 | 0.03 | 0.03 | 0.17 | 0.18 | 0.06 | 0.03 |
| Lincoln-Hitler | 0.29 | 0.32 | 0.11 | 0.11 | 0.08 | 0.11 | 0.11 | 0.11 |
| Personality - Agreeableness | 0.31 | 0.31 | 0.13 | 0.21 | 0.05 | 0.13 | 0.15 | 0.15 |
| Personality - Conscientiousness | 0.38 | 0.41 | 0.03 | 0.00 | 0.10 | 0.10 | 0.10 | 0.21 |
| Personality - Extraversion | 0.15 | 0.15 | 0.09 | 0.06 | 0.21 | 0.21 | 0.03 | 0.00 |
| Personality - Neuroticism | 0.25 | 0.28 | 0.06 | 0.03 | 0.22 | 0.16 | 0.16 | 0.16 |
| Personality - Openness | 0.12 | 0.09 | 0.12 | 0.06 | 0.06 | 0.03 | 0.15 | 0.12 |
| Race (1) | 0.16 | 0.18 | 0.02 | 0.02 | 0.16 | 0.16 | 0.09 | 0.09 |
| Race (2) | 0.36 | 0.43 | 0.07 | 0.11 | 0.14 | 0.20 | 0.11 | 0.09 |
| Religion | 0.17 | 0.17 | 0.10 | 0.07 | 0.17 | 0.20 | 0.13 | 0.07 |
| Rich-Poor | 0.18 | 0.25 | 0.11 | 0.08 | 0.12 | 0.19 | 0.07 | 0.11 |
| Sexuality (1) | 0.42 | 0.50 | 0.04 | 0.15 | 0.27 | 0.35 | 0.15 | 0.12 |
| Sexuality (2) | 0.58 | 0.58 | 0.37 | 0.37 | 0.42 | 0.47 | 0.26 | 0.37 |
| Shapes & colors (1) | 0.18 | 0.27 | 0.18 | 0.18 | 0.00 | 0.09 | 0.00 | 0.00 |
| Shapes & colors (2) | 0.18 | 0.36 | 0.09 | 0.09 | 0.00 | 0.00 | 0.09 | 0.18 |
| Shapes & colors (3) | 0.10 | 0.10 | 0.10 | 0.10 | 0.30 | 0.30 | 0.20 | 0.30 |
| Shapes & colors (4) | 0.25 | 0.38 | 0.00 | 0.00 | 0.00 | 0.12 | 0.25 | 0.12 |
| Shapes & colors (5) | 0.36 | 0.50 | 0.00 | 0.05 | 0.00 | 0.09 | 0.00 | 0.09 |
| Shapes & colors (6) | 0.23 | 0.23 | 0.12 | 0.19 | 0.08 | 0.08 | 0.27 | 0.31 |
| Shapes & colors (7) | 0.57 | 0.50 | 0.14 | 0.21 | 0.29 | 0.29 | 0.36 | 0.36 |
| Stigma - ADHD | 0.16 | 0.27 | 0.08 | 0.08 | 0.13 | 0.16 | 0.11 | 0.18 |
| Stigma - PTSD | 0.30 | 0.31 | 0.11 | 0.07 | 0.09 | 0.20 | 0.11 | 0.09 |
| Stigma - Schizophrenia | 0.13 | 0.22 | 0.11 | 0.11 | 0.22 | 0.30 | 0.15 | 0.13 |
| Valenced words | 0.41 | 0.49 | 0.14 | 0.11 | 0.08 | 0.14 | 0.08 | 0.11 |
p_diff_zero_scoring <-
data_diff_zero %>%
filter(method == "bca") %>%
mutate(domain = fct_reorder(domain, proportion_diff_zero, .fun = mean)) %>%
ggplot(aes(proportion_diff_zero, domain, color = DV_type, shape = DV_type)) +
geom_linerangeh(aes(xmin = proportion_diff_zero - sqrt(variance)*1.96,
xmax = proportion_diff_zero + sqrt(variance)*1.96),
position = position_dodge(width = 1)) +
geom_point(position = position_dodge(width = 1)) +
scale_shape_manual(labels = c("*D* scores", "PI scores"), values = c(15, 16)) +
scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("*D* scores", "PI scores")) +
mdthemes::md_theme_linedraw() +
facet_wrap(~ trial_type, ncol = 4) +
labs(x = "Proportion of scores different from zero point",
y = "",
color = "Scoring method",
shape = "Scoring method") +
theme(legend.position = "top")
p_diff_zero_scoringp_diff_zero_method <-
data_diff_zero %>%
filter(DV_type == "*D* scores") %>%
mutate(domain = fct_reorder(domain, proportion_diff_zero, .fun = mean)) %>%
ggplot(aes(proportion_diff_zero, domain, color = method, shape = method)) +
geom_linerangeh(aes(xmin = proportion_diff_zero - sqrt(variance)*1.96,
xmax = proportion_diff_zero + sqrt(variance)*1.96),
position = position_dodge(width = 1)) +
geom_point(position = position_dodge(width = 1)) +
scale_shape_manual(labels = c("BCA", "Basic", "Normal", "Percentile"), values = c(15, 16, 17, 7)) +
scale_color_viridis_d(begin = 0.2, end = 0.8, labels = c("BCA", "Basic", "Normal", "Percentile")) +
mdthemes::md_theme_linedraw() +
facet_wrap(~ trial_type, ncol = 4) +
labs(x = "Proportion of scores different from zero point",
y = "",
color = "Bootstrapping method",
shape = "Bootstrapping method") +
theme(legend.position = "top")
p_diff_zero_method# fit model
fit_diff_zero_scoring <-
data_diff_zero %>%
filter(method == "bca") %>%
lmer(proportion_diff_zero ~ 1 + trial_type * DV_type + (trial_type | domain),
weights = 1/variance,
data = .)
# results table
tab_model(fit_diff_zero_scoring,
string.p = "p (corrected)",
ci.hyphen = ", ",
#emph.p = FALSE,
p.adjust = "holm")| proportion diff zero | |||
|---|---|---|---|
| Predictors | Estimates | CI | p (corrected) |
| (Intercept) | 0.24 | 0.19, 0.28 | <0.001 |
| trial_type [Trial type 2] | -0.14 | -0.18, -0.09 | <0.001 |
| trial_type [Trial type 3] | -0.10 | -0.15, -0.05 | <0.001 |
| trial_type [Trial type 4] | -0.11 | -0.15, -0.06 | <0.001 |
| DV_type PI scores | 0.04 | 0.02, 0.06 | 0.001 |
|
trial_type [Trial type 2] * DV_type PI scores |
-0.04 | -0.06, -0.02 | 0.001 |
|
trial_type [Trial type 3] * DV_type PI scores |
-0.04 | -0.06, -0.02 | 0.001 |
|
trial_type [Trial type 4] * DV_type PI scores |
-0.04 | -0.06, -0.02 | 0.001 |
| Random Effects | |||
| σ2 | 0.42 | ||
| τ00 domain | 0.01 | ||
| τ11 domain.trial_typeTrial type 2 | 0.01 | ||
| τ11 domain.trial_typeTrial type 3 | 0.02 | ||
| τ11 domain.trial_typeTrial type 4 | 0.02 | ||
| ρ01 | -0.85 | ||
| -0.76 | |||
| -0.76 | |||
| ICC | 0.02 | ||
| N domain | 35 | ||
| Observations | 280 | ||
| Marginal R2 / Conditional R2 | 0.008 / 0.029 | ||
# summary
summary(fit_diff_zero_scoring)## Linear mixed model fit by REML ['lmerMod']
## Formula: proportion_diff_zero ~ 1 + trial_type * DV_type + (trial_type |
## domain)
## Data: .
## Weights: 1/variance
##
## REML criterion at convergence: -669.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.45472 -0.29090 0.08845 0.66833 2.23314
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## domain (Intercept) 0.01478 0.1216
## trial_typeTrial type 2 0.01408 0.1187 -0.85
## trial_typeTrial type 3 0.01864 0.1365 -0.76 0.74
## trial_typeTrial type 4 0.01698 0.1303 -0.76 0.78 0.73
## Residual 0.42016 0.6482
## Number of obs: 280, groups: domain, 35
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.235973 0.021956 10.747
## trial_typeTrial type 2 -0.136759 0.021897 -6.246
## trial_typeTrial type 3 -0.101478 0.024771 -4.097
## trial_typeTrial type 4 -0.107244 0.023846 -4.497
## DV_typePI scores 0.036650 0.009736 3.764
## trial_typeTrial type 2:DV_typePI scores -0.036574 0.009778 -3.740
## trial_typeTrial type 3:DV_typePI scores -0.036282 0.009779 -3.710
## trial_typeTrial type 4:DV_typePI scores -0.036576 0.009778 -3.740
##
## Correlation of Fixed Effects:
## (Intr) tr_Tt2 tr_Tt3 tr_Tt4 DV_PIs t_Tt2s t_Tt3s
## trl_typTrt2 -0.854
## trl_typTrt3 -0.770 0.743
## trl_typTrt4 -0.770 0.775 0.731
## DV_typPIscr -0.210 0.211 0.186 0.194
## t_Tt2:DV_Ps 0.210 -0.212 -0.186 -0.193 -0.996
## t_Tt3:DV_Ps 0.210 -0.210 -0.187 -0.193 -0.996 0.991
## t_Tt4:DV_Ps 0.210 -0.210 -0.186 -0.195 -0.996 0.991 0.991
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00279646 (tol = 0.002, component 1)
# plot RE effects
plot_model(fit_diff_zero_scoring,
colors = "bw",
type = "re") +
mdthemes::md_theme_linedraw() +
ggtitle("")# extract re Tau
results_re_tau_diff_zero_scoring <- fit_diff_zero_scoring %>%
merTools::REsdExtract() %>%
as_tibble(rownames = "trial_type") %>%
rename(tau = value) %>%
mutate(trial_type = case_when(trial_type == "domain.(Intercept)" ~ "Trial type 1",
trial_type == "domain.trial_typeTrial type 2" ~ "Trial type 2",
trial_type == "domain.trial_typeTrial type 3" ~ "Trial type 3",
trial_type == "domain.trial_typeTrial type 4" ~ "Trial type 4"))
# extract marginal means
results_emm_diff_zero_scoring <-
summary(emmeans(fit_diff_zero_scoring, ~ DV_type | trial_type)) %>%
dplyr::select(DV_type, trial_type, estimate = emmean, se = SE, ci_lower = lower.CL, ci_upper = upper.CL) %>%
mutate(trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4"))
# combine
results_diff_zero_scoring <- results_emm_diff_zero_scoring %>%
full_join(results_re_tau_diff_zero_scoring, by = "trial_type") %>%
mutate(cr_lower = estimate - (1.96 * sqrt(se^2 + tau^2)), # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
cr_upper = estimate + (1.96 * sqrt(se^2 + tau^2)))
# save to disk
write_rds(results_diff_zero_scoring, "models/results_diff_zero_scoring.rds")
# tests
emmeans(fit_diff_zero_scoring, list(pairwise ~ DV_type | trial_type), adjust = "bonferroni")## $`emmeans of DV_type | trial_type`
## trial_type = Trial type 1:
## DV_type emmean SE df lower.CL upper.CL
## *D* scores 0.2360 0.0222 6625 0.1863 0.286
## PI scores 0.2726 0.0223 6628 0.2227 0.323
##
## trial_type = Trial type 2:
## DV_type emmean SE df lower.CL upper.CL
## *D* scores 0.0992 0.0119 25313 0.0725 0.126
## PI scores 0.0993 0.0119 25318 0.0726 0.126
##
## trial_type = Trial type 3:
## DV_type emmean SE df lower.CL upper.CL
## *D* scores 0.1345 0.0162 13980 0.0982 0.171
## PI scores 0.1349 0.0162 13982 0.0986 0.171
##
## trial_type = Trial type 4:
## DV_type emmean SE df lower.CL upper.CL
## *D* scores 0.1287 0.0157 13415 0.0935 0.164
## PI scores 0.1288 0.0157 13416 0.0935 0.164
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
##
## $`pairwise differences of DV_type | trial_type`
## trial_type = Trial type 1:
## contrast estimate SE df t.ratio p.value
## *D* scores - PI scores -0.0366504 0.009737 8979469 -3.764 0.0002
##
## trial_type = Trial type 2:
## contrast estimate SE df t.ratio p.value
## *D* scores - PI scores -0.0000768 0.000906 119717151791 -0.085 0.9325
##
## trial_type = Trial type 3:
## contrast estimate SE df t.ratio p.value
## *D* scores - PI scores -0.0003681 0.000910 117400475158 -0.404 0.6859
##
## trial_type = Trial type 4:
## contrast estimate SE df t.ratio p.value
## *D* scores - PI scores -0.0000744 0.000906 119551896105 -0.082 0.9346
##
## Degrees-of-freedom method: kenward-roger
emmeans(fit_diff_zero_scoring, list(pairwise ~ trial_type | DV_type), adjust = "bonferroni")## $`emmeans of trial_type | DV_type`
## DV_type = *D* scores:
## trial_type emmean SE df lower.CL upper.CL
## Trial type 1 0.2360 0.0222 6625 0.1806 0.291
## Trial type 2 0.0992 0.0119 25313 0.0695 0.129
## Trial type 3 0.1345 0.0162 13980 0.0941 0.175
## Trial type 4 0.1287 0.0157 13415 0.0894 0.168
##
## DV_type = PI scores:
## trial_type emmean SE df lower.CL upper.CL
## Trial type 1 0.2726 0.0223 6628 0.2170 0.328
## Trial type 2 0.0993 0.0119 25318 0.0695 0.129
## Trial type 3 0.1349 0.0162 13982 0.0945 0.175
## Trial type 4 0.1288 0.0157 13416 0.0895 0.168
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
##
## $`pairwise differences of trial_type | DV_type`
## DV_type = *D* scores:
## contrast estimate SE df t.ratio p.value
## Trial type 1 - Trial type 2 0.13676 0.0220 13560 6.203 <.0001
## Trial type 1 - Trial type 3 0.10148 0.0249 13817 4.069 0.0003
## Trial type 1 - Trial type 4 0.10724 0.0240 15081 4.468 <.0001
## Trial type 2 - Trial type 3 -0.03528 0.0170 46982 -2.074 0.2283
## Trial type 2 - Trial type 4 -0.02951 0.0155 61762 -1.904 0.3416
## Trial type 3 - Trial type 4 0.00577 0.0179 50978 0.322 1.0000
##
## DV_type = PI scores:
## contrast estimate SE df t.ratio p.value
## Trial type 1 - Trial type 2 0.17333 0.0222 13508 7.822 <.0001
## Trial type 1 - Trial type 3 0.13776 0.0250 13768 5.501 <.0001
## Trial type 1 - Trial type 4 0.14382 0.0241 15020 5.966 <.0001
## Trial type 2 - Trial type 3 -0.03557 0.0170 46980 -2.091 0.2192
## Trial type 2 - Trial type 4 -0.02951 0.0155 61768 -1.903 0.3419
## Trial type 3 - Trial type 4 0.00606 0.0179 50977 0.338 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: bonferroni method for 6 tests
# plot
dodge_width <- 0.8
p_prop_nonzero_scoring <-
ggplot(results_diff_zero_scoring, aes(fct_rev(trial_type), estimate, color = DV_type, shape = DV_type, group = DV_type)) +
geom_linerange(aes(ymin = cr_lower, ymax = cr_upper), position = position_dodge(width = dodge_width), linetype = "dotted") +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1, position = position_dodge(width = dodge_width)) +
# geom_line(position = position_dodge(width = dodge_width)) +
geom_point(position = position_dodge(width = dodge_width), size = 2.5) +
mdthemes::md_theme_linedraw() +
scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Worse)", "0.25", "0.50", "0.75", "1.00<br/>(Better)")) +
scale_color_viridis_d(alpha = 1, begin = 0.3, end = 0.7, labels = c("*D* scores", "PI scores")) +
labs(shape = "Scoring method",
color = "Scoring method",
x = "",
y = "Proportion of participants with non-zero scores<br/>") +
theme(legend.position = "top") +
coord_flip(ylim = c(0, 1))
p_prop_nonzero_scoring# fit model
fit_diff_zero_method <-
data_diff_zero %>%
filter(DV_type == "*D* scores") %>%
lmer(proportion_diff_zero ~ 1 + trial_type * method + (trial_type | domain),
weights = 1/variance,
data = .)
# results table
tab_model(fit_diff_zero_method,
string.p = "p (corrected)",
ci.hyphen = ", ",
#emph.p = FALSE,
p.adjust = "holm")| proportion diff zero | |||
|---|---|---|---|
| Predictors | Estimates | CI | p (corrected) |
| (Intercept) | 0.31 | 0.26, 0.35 | <0.001 |
| trial_type [Trial type 2] | -0.19 | -0.24, -0.15 | <0.001 |
| trial_type [Trial type 3] | -0.16 | -0.22, -0.10 | <0.001 |
| trial_type [Trial type 4] | -0.13 | -0.18, -0.08 | <0.001 |
| method [normal] | 0.03 | 0.01, 0.05 | 0.007 |
| method [basic] | 0.06 | 0.05, 0.08 | <0.001 |
| method [bca] | -0.07 | -0.08, -0.05 | <0.001 |
|
trial_type [Trial type 2] * method [normal] |
0.03 | 0.00, 0.05 | 0.068 |
|
trial_type [Trial type 3] * method [normal] |
0.04 | 0.01, 0.06 | 0.010 |
|
trial_type [Trial type 4] * method [normal] |
-0.01 | -0.03, 0.02 | 0.996 |
|
trial_type [Trial type 2] * method [basic] |
0.03 | 0.01, 0.06 | 0.024 |
|
trial_type [Trial type 3] * method [basic] |
0.04 | 0.02, 0.06 | 0.009 |
|
trial_type [Trial type 4] * method [basic] |
-0.01 | -0.03, 0.02 | 0.996 |
|
trial_type [Trial type 2] * method [bca] |
0.07 | 0.05, 0.08 | <0.001 |
|
trial_type [Trial type 3] * method [bca] |
0.07 | 0.05, 0.08 | <0.001 |
|
trial_type [Trial type 4] * method [bca] |
0.02 | -0.00, 0.04 | 0.151 |
| Random Effects | |||
| σ2 | 0.33 | ||
| τ00 domain | 0.02 | ||
| τ11 domain.trial_typeTrial type 2 | 0.02 | ||
| τ11 domain.trial_typeTrial type 3 | 0.03 | ||
| τ11 domain.trial_typeTrial type 4 | 0.02 | ||
| ρ01 | -0.84 | ||
| -0.78 | |||
| -0.71 | |||
| ICC | 0.03 | ||
| N domain | 35 | ||
| Observations | 560 | ||
| Marginal R2 / Conditional R2 | 0.016 / 0.048 | ||
# summary
summary(fit_diff_zero_method)## Linear mixed model fit by REML ['lmerMod']
## Formula: proportion_diff_zero ~ 1 + trial_type * method + (trial_type |
## domain)
## Data: .
## Weights: 1/variance
##
## REML criterion at convergence: -1517.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2582 -0.3883 0.1004 0.5687 3.1523
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## domain (Intercept) 0.01746 0.1321
## trial_typeTrial type 2 0.01680 0.1296 -0.84
## trial_typeTrial type 3 0.02650 0.1628 -0.78 0.80
## trial_typeTrial type 4 0.02095 0.1447 -0.71 0.71 0.73
## Residual 0.33022 0.5746
## Number of obs: 560, groups: domain, 35
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.308929 0.023434 13.183
## trial_typeTrial type 2 -0.192111 0.023318 -8.239
## trial_typeTrial type 3 -0.160350 0.028678 -5.591
## trial_typeTrial type 4 -0.127595 0.026081 -4.892
## methodnormal 0.031708 0.009494 3.340
## methodbasic 0.064431 0.009627 6.693
## methodbca -0.067421 0.008871 -7.600
## trial_typeTrial type 2:methodnormal 0.027250 0.011405 2.389
## trial_typeTrial type 3:methodnormal 0.036709 0.011693 3.140
## trial_typeTrial type 4:methodnormal -0.008180 0.012075 -0.677
## trial_typeTrial type 2:methodbasic 0.033191 0.011771 2.820
## trial_typeTrial type 3:methodbasic 0.038777 0.012038 3.221
## trial_typeTrial type 4:methodbasic -0.005323 0.012405 -0.429
## trial_typeTrial type 2:methodbca 0.066640 0.008908 7.481
## trial_typeTrial type 3:methodbca 0.066750 0.008908 7.493
## trial_typeTrial type 4:methodbca 0.021179 0.010817 1.958
# plot re
plot_model(fit_diff_zero_method,
colors = "bw",
type = "re") +
mdthemes::md_theme_linedraw() +
ggtitle("")# extract re Tau
results_re_tau_diff_zero_method <- fit_diff_zero_method %>%
merTools::REsdExtract() %>%
as_tibble(rownames = "trial_type") %>%
rename(tau = value) %>%
mutate(trial_type = case_when(trial_type == "domain.(Intercept)" ~ "Trial type 1",
trial_type == "domain.trial_typeTrial type 2" ~ "Trial type 2",
trial_type == "domain.trial_typeTrial type 3" ~ "Trial type 3",
trial_type == "domain.trial_typeTrial type 4" ~ "Trial type 4"))
# extract marginal means
results_emm_diff_zero_method <-
summary(emmeans(fit_diff_zero_method, ~ method | trial_type)) %>%
dplyr::select(method, trial_type, estimate = emmean, se = SE, ci_lower = lower.CL, ci_upper = upper.CL) %>%
mutate(trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4"))
# combine
results_diff_zero_method <- results_emm_diff_zero_method %>%
full_join(results_re_tau_diff_zero_method, by = "trial_type") %>%
mutate(cr_lower = estimate - (1.96 * sqrt(se^2 + tau^2)), # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
cr_upper = estimate + (1.96 * sqrt(se^2 + tau^2)))
# save to disk
write_rds(results_diff_zero_method, "models/results_diff_zero_method.rds")
# tests
emmeans(fit_diff_zero_method, list(pairwise ~ method | trial_type), adjust = "bonferroni")## $`emmeans of method | trial_type`
## trial_type = Trial type 1:
## method emmean SE df lower.CL upper.CL
## percent 0.309 0.0240 1080 0.2490 0.369
## normal 0.341 0.0240 1086 0.2806 0.401
## basic 0.373 0.0241 1091 0.3132 0.434
## bca 0.242 0.0238 1063 0.1820 0.301
##
## trial_type = Trial type 2:
## method emmean SE df lower.CL upper.CL
## percent 0.117 0.0132 3220 0.0838 0.150
## normal 0.176 0.0140 3885 0.1408 0.211
## basic 0.214 0.0142 4040 0.1790 0.250
## bca 0.116 0.0132 3215 0.0830 0.149
##
## trial_type = Trial type 3:
## method emmean SE df lower.CL upper.CL
## percent 0.149 0.0179 1946 0.1039 0.193
## normal 0.217 0.0186 2180 0.1706 0.263
## basic 0.252 0.0187 2220 0.2050 0.299
## bca 0.148 0.0179 1945 0.1032 0.193
##
## trial_type = Trial type 4:
## method emmean SE df lower.CL upper.CL
## percent 0.181 0.0191 1673 0.1335 0.229
## normal 0.205 0.0193 1704 0.1567 0.253
## basic 0.240 0.0194 1735 0.1920 0.289
## bca 0.135 0.0189 1619 0.0879 0.182
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
##
## $`pairwise differences of method | trial_type`
## trial_type = Trial type 1:
## contrast estimate SE df t.ratio p.value
## percent - normal -0.031708 0.009494 18390245 -3.340 0.0050
## percent - basic -0.064431 0.009627 17350004 -6.693 <.0001
## percent - bca 0.067421 0.008872 23853769 7.600 <.0001
## normal - basic -0.032723 0.009766 16426027 -3.351 0.0048
## normal - bca 0.099129 0.009032 22034114 10.975 <.0001
## basic - bca 0.131852 0.009176 20512833 14.369 <.0001
##
## trial_type = Trial type 2:
## contrast estimate SE df t.ratio p.value
## percent - normal -0.058958 0.006320 52184590 -9.328 <.0001
## percent - basic -0.097622 0.006774 38249409 -14.411 <.0001
## percent - bca 0.000781 0.000806 330766001965 0.969 1.0000
## normal - basic -0.038664 0.008135 33864472 -4.753 <.0001
## normal - bca 0.059738 0.006308 52119682 9.471 <.0001
## basic - bca 0.098403 0.006763 38176259 14.551 <.0001
##
## trial_type = Trial type 3:
## contrast estimate SE df t.ratio p.value
## percent - normal -0.068417 0.006826 38328005 -10.023 <.0001
## percent - basic -0.103208 0.007228 29755644 -14.279 <.0001
## percent - bca 0.000671 0.000807 335130816874 0.831 1.0000
## normal - basic -0.034791 0.008684 26144671 -4.006 0.0004
## normal - bca 0.069088 0.006816 38292899 10.136 <.0001
## basic - bca 0.103879 0.007219 29717292 14.390 <.0001
##
## trial_type = Trial type 4:
## contrast estimate SE df t.ratio p.value
## percent - normal -0.023528 0.007462 47629729 -3.153 0.0097
## percent - basic -0.059108 0.007823 38385761 -7.555 <.0001
## percent - bca 0.046242 0.006189 94703959 7.472 <.0001
## normal - basic -0.035580 0.008115 34205153 -4.384 0.0001
## normal - bca 0.069770 0.006756 62785077 10.326 <.0001
## basic - bca 0.105350 0.007200 45936913 14.631 <.0001
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: bonferroni method for 6 tests
emmeans(fit_diff_zero_method, list(pairwise ~ trial_type | method), adjust = "bonferroni")## $`emmeans of trial_type | method`
## method = percent:
## trial_type emmean SE df lower.CL upper.CL
## Trial type 1 0.309 0.0240 1080 0.2490 0.369
## Trial type 2 0.117 0.0132 3220 0.0838 0.150
## Trial type 3 0.149 0.0179 1946 0.1039 0.193
## Trial type 4 0.181 0.0191 1673 0.1335 0.229
##
## method = normal:
## trial_type emmean SE df lower.CL upper.CL
## Trial type 1 0.341 0.0240 1086 0.2806 0.401
## Trial type 2 0.176 0.0140 3885 0.1408 0.211
## Trial type 3 0.217 0.0186 2180 0.1706 0.263
## Trial type 4 0.205 0.0193 1704 0.1567 0.253
##
## method = basic:
## trial_type emmean SE df lower.CL upper.CL
## Trial type 1 0.373 0.0241 1091 0.3132 0.434
## Trial type 2 0.214 0.0142 4040 0.1790 0.250
## Trial type 3 0.252 0.0187 2220 0.2050 0.299
## Trial type 4 0.240 0.0194 1735 0.1920 0.289
##
## method = bca:
## trial_type emmean SE df lower.CL upper.CL
## Trial type 1 0.242 0.0238 1063 0.1820 0.301
## Trial type 2 0.116 0.0132 3215 0.0830 0.149
## Trial type 3 0.148 0.0179 1945 0.1032 0.193
## Trial type 4 0.135 0.0189 1619 0.0879 0.182
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
##
## $`pairwise differences of trial_type | method`
## method = percent:
## contrast estimate SE df t.ratio p.value
## Trial type 1 - Trial type 2 0.1921 0.0237 2123 8.109 <.0001
## Trial type 1 - Trial type 3 0.1604 0.0292 1674 5.493 <.0001
## Trial type 1 - Trial type 4 0.1276 0.0265 2380 4.819 <.0001
## Trial type 2 - Trial type 3 -0.0318 0.0175 6443 -1.814 0.4187
## Trial type 2 - Trial type 4 -0.0645 0.0191 6029 -3.371 0.0045
## Trial type 3 - Trial type 4 -0.0328 0.0208 5439 -1.578 0.6882
##
## method = normal:
## contrast estimate SE df t.ratio p.value
## Trial type 1 - Trial type 2 0.1649 0.0242 2284 6.819 <.0001
## Trial type 1 - Trial type 3 0.1236 0.0297 1756 4.170 0.0002
## Trial type 1 - Trial type 4 0.1358 0.0266 2413 5.102 <.0001
## Trial type 2 - Trial type 3 -0.0412 0.0188 7947 -2.197 0.1683
## Trial type 2 - Trial type 4 -0.0291 0.0198 6724 -1.469 0.8512
## Trial type 3 - Trial type 4 0.0121 0.0215 5975 0.566 1.0000
##
## method = basic:
## contrast estimate SE df t.ratio p.value
## Trial type 1 - Trial type 2 0.1589 0.0243 2331 6.529 <.0001
## Trial type 1 - Trial type 3 0.1216 0.0298 1774 4.082 0.0003
## Trial type 1 - Trial type 4 0.1329 0.0267 2445 4.969 <.0001
## Trial type 2 - Trial type 3 -0.0373 0.0190 8242 -1.961 0.2996
## Trial type 2 - Trial type 4 -0.0260 0.0201 6951 -1.296 1.0000
## Trial type 3 - Trial type 4 0.0113 0.0217 6124 0.523 1.0000
##
## method = bca:
## contrast estimate SE df t.ratio p.value
## Trial type 1 - Trial type 2 0.1255 0.0235 2094 5.337 <.0001
## Trial type 1 - Trial type 3 0.0936 0.0290 1660 3.223 0.0078
## Trial type 1 - Trial type 4 0.1064 0.0261 2312 4.073 0.0003
## Trial type 2 - Trial type 3 -0.0319 0.0175 6435 -1.821 0.4121
## Trial type 2 - Trial type 4 -0.0191 0.0189 5901 -1.010 1.0000
## Trial type 3 - Trial type 4 0.0128 0.0205 5314 0.625 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: bonferroni method for 6 tests
# plot
p_prop_nonzero_method <-
ggplot(results_diff_zero_method, aes(fct_rev(trial_type), estimate, color = method, shape = method, group = method)) +
geom_linerange(aes(ymin = cr_lower, ymax = cr_upper), position = position_dodge(width = dodge_width), linetype = "dotted") +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1, position = position_dodge(width = dodge_width)) +
# geom_line(position = position_dodge(width = dodge_width)) +
geom_point(position = position_dodge(width = dodge_width), size = 2.5) +
mdthemes::md_theme_linedraw() +
scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Worse)", "0.25", "0.50", "0.75", "1.00<br/>(Better)")) +
scale_shape_discrete(labels = c("BCA", "Basic", "Normal", "Percentile")) +
scale_color_viridis_d(alpha = 1, begin = 0.3, end = 0.7, labels = c("BCA", "Basic", "Normal", "Percentile")) +
labs(shape = "Bootstrapping method",
color = "Bootstrapping method",
x = "",
y = "Proportion of participants with non-zero scores") +
theme(legend.position = "top") +
coord_flip(ylim = c(0, 1))
p_prop_nonzero_methodWithin domain and trial type
Many have argued that the zero point is arbitrary and not a useful reference point. Instead of asking “what proportion of D/PI scores are different from zero?”, we could also ask “what proportion of D/PI scores are different from one another?”
# helper function to apply workflow to each resample
discriminability <- function(data, i) {
data_with_indexes <- data[i,] # boot function requires data and index
estimate <- data_with_indexes$estimate
ci_lower <- data_with_indexes$ci_lower
ci_upper <- data_with_indexes$ci_upper
n_estimate <- length(estimate)
n_ci_lower <- length(ci_lower)
n_ci_upper <- length(ci_upper)
r_estimate <- sum(rank(c(estimate, ci_lower))[1:n_estimate])
r_ci_upper <- sum(rank(c(ci_upper, estimate))[1:n_ci_upper])
prob_estimate_inferior_to_ci_lower <- 1 - (r_estimate / n_estimate - (n_estimate + 1) / 2) / n_ci_lower
prob_estimate_superior_to_ci_upper <- 1 - (r_ci_upper / n_ci_upper - (n_ci_upper + 1) / 2) / n_estimate
probability_estimates_outside_cis <- (prob_estimate_inferior_to_ci_lower + prob_estimate_superior_to_ci_upper)
return(probability_estimates_outside_cis)
}
bootstrap_discriminability <- function(data){
require(dplyr)
require(boot)
fit <-
boot::boot(data = data,
statistic = discriminability,
R = n_boots,
sim = "ordinary",
stype = "i",
parallel = "multicore",
ncpus = parallel::detectCores())
results <- boot::boot.ci(fit, conf = 0.95, type = c("norm","basic", "perc", "bca"))
output <-
tibble(method = c("normal", "basic", "percent", "bca"),
estimate = rep(fit$t0, 4),
ci_lower = c(results$normal[2], results$basic[4], results$percent[4], results$bca[4]),
ci_upper = c(results$normal[3], results$basic[5], results$percent[5], results$bca[5]))
return(output)
}# bootstrapping has a long execution time, so load saved values if they've already been calculated
if(file.exists("models/data_discriminability_D.rds")) {
data_discriminability_D <- read_rds("models/data_discriminability_D.rds")
} else {
# bootstrap D scores
data_discriminability_D <- data_estimates_D %>%
select(unique_id, domain, trial_type, estimate, ci_upper, ci_lower) %>%
group_by(domain, trial_type) %>%
do(bootstrap_discriminability(data = .)) %>%
ungroup() %>%
mutate(proportion_discriminable = estimate,
variance = ((ci_upper - ci_lower)/(1.96*2))^2) %>%
mutate(domain = as.factor(domain),
method = fct_relevel(method, "bca", "basic", "normal", "percent"),
trial_type = fct_relevel(trial_type, "tt1", "tt2", "tt3", "tt4"),
DV_type = "*D* scores")
# save to disk
write_rds(data_discriminability_D, "models/data_discriminability_D.rds", compress = "gz")
}# bootstrapping has a long execution time, so load saved values if they've already been calculated
if(file.exists("models/data_discriminability_PI.rds")) {
data_discriminability_PI <- read_rds("models/data_discriminability_PI.rds")
} else {
# bootstrap D scores
data_discriminability_PI <- data_estimates_PI %>%
select(unique_id, domain, trial_type, estimate, ci_upper, ci_lower) %>%
group_by(domain, trial_type) %>%
do(bootstrap_discriminability(data = .)) %>%
ungroup() %>%
mutate(proportion_discriminable = estimate,
variance = ((ci_upper - ci_lower)/(1.96*2))^2) %>%
mutate(domain = as.factor(domain),
method = fct_relevel(method, "bca", "basic", "normal", "percent"),
trial_type = fct_relevel(trial_type, "tt1", "tt2", "tt3", "tt4"),
DV_type = "PI scores")
# save to disk
write_rds(data_discriminability_PI, "models/data_discriminability_PI.rds", compress = "gz")
}# combine
data_discriminability_combined_scoring <-
bind_rows(data_discriminability_D,
data_discriminability_PI) %>%
filter(method == "bca") %>%
mutate(trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
trial_type == "tt2" ~ "Trial type 2",
trial_type == "tt3" ~ "Trial type 3",
trial_type == "tt4" ~ "Trial type 4"))
# fit meta analytic model
fit_disciminability_scoring <-
lmer(proportion_discriminable ~ 1 + trial_type * DV_type + (trial_type | domain),
weights = 1/variance,
data = data_discriminability_combined_scoring)
# results table
tab_model(fit_disciminability_scoring,
string.p = "p (corrected)",
ci.hyphen = ", ",
#emph.p = FALSE,
p.adjust = "holm")| proportion discriminable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p (corrected) |
| (Intercept) | 0.26 | 0.23, 0.29 | <0.001 |
| trial_type [Trial type 2] | -0.02 | -0.05, 0.01 | 1.000 |
| trial_type [Trial type 3] | 0.01 | -0.03, 0.04 | 1.000 |
| trial_type [Trial type 4] | -0.00 | -0.03, 0.03 | 1.000 |
| DV_type PI scores | -0.00 | -0.01, 0.01 | 1.000 |
|
trial_type [Trial type 2] * DV_type PI scores |
-0.01 | -0.02, 0.01 | 1.000 |
|
trial_type [Trial type 3] * DV_type PI scores |
-0.01 | -0.02, 0.01 | 1.000 |
|
trial_type [Trial type 4] * DV_type PI scores |
-0.00 | -0.02, 0.01 | 1.000 |
| Random Effects | |||
| σ2 | 0.94 | ||
| τ00 domain | 0.01 | ||
| τ11 domain.trial_typeTrial type 2 | 0.01 | ||
| τ11 domain.trial_typeTrial type 3 | 0.01 | ||
| τ11 domain.trial_typeTrial type 4 | 0.01 | ||
| ρ01 | -0.62 | ||
| -0.42 | |||
| -0.65 | |||
| ICC | 0.01 | ||
| N domain | 35 | ||
| Observations | 280 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.007 | ||
# plot RE effects
plot_model(fit_disciminability_scoring,
colors = "bw",
type = "re") +
mdthemes::md_theme_linedraw() +
ggtitle("")# extract re Tau
results_re_tau_disciminability_scoring <- fit_disciminability_scoring %>%
merTools::REsdExtract() %>%
as_tibble(rownames = "trial_type") %>%
rename(tau = value) %>%
mutate(trial_type = case_when(trial_type == "domain.(Intercept)" ~ "Trial type 1",
trial_type == "domain.trial_typeTrial type 2" ~ "Trial type 2",
trial_type == "domain.trial_typeTrial type 3" ~ "Trial type 3",
trial_type == "domain.trial_typeTrial type 4" ~ "Trial type 4"))
# extract marginal means
results_emm_disciminability_scoring <-
summary(emmeans(fit_disciminability_scoring, ~ DV_type | trial_type)) %>%
dplyr::select(DV_type, trial_type, estimate = emmean, se = SE, ci_lower = lower.CL, ci_upper = upper.CL) %>%
mutate(trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4"))
# combine
results_disciminability_scoring <- results_emm_disciminability_scoring %>%
full_join(results_re_tau_disciminability_scoring, by = "trial_type") %>%
mutate(cr_lower = estimate - (1.96 * sqrt(se^2 + tau^2)), # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
cr_upper = estimate + (1.96 * sqrt(se^2 + tau^2)))
# save to disk
write_rds(results_disciminability_scoring, "models/results_disciminability_scoring.rds")
# tests
emmeans(fit_disciminability_scoring, list(pairwise ~ DV_type | trial_type), adjust = "bonferroni")## $`emmeans of DV_type | trial_type`
## trial_type = Trial type 1:
## DV_type emmean SE df lower.CL upper.CL
## *D* scores 0.257 0.0147 84755 0.224 0.290
## PI scores 0.257 0.0148 85049 0.224 0.290
##
## trial_type = Trial type 2:
## DV_type emmean SE df lower.CL upper.CL
## *D* scores 0.236 0.0132 72585 0.207 0.266
## PI scores 0.230 0.0132 72895 0.200 0.259
##
## trial_type = Trial type 3:
## DV_type emmean SE df lower.CL upper.CL
## *D* scores 0.267 0.0174 41481 0.227 0.306
## PI scores 0.261 0.0175 41607 0.222 0.300
##
## trial_type = Trial type 4:
## DV_type emmean SE df lower.CL upper.CL
## *D* scores 0.257 0.0129 104739 0.228 0.286
## PI scores 0.254 0.0129 104745 0.225 0.283
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
##
## $`pairwise differences of DV_type | trial_type`
## trial_type = Trial type 1:
## contrast estimate SE df t.ratio p.value
## *D* scores - PI scores 0.000148 0.00606 302334465 0.024 0.9805
##
## trial_type = Trial type 2:
## contrast estimate SE df t.ratio p.value
## *D* scores - PI scores 0.006683 0.00579 361738352 1.154 0.2486
##
## trial_type = Trial type 3:
## contrast estimate SE df t.ratio p.value
## *D* scores - PI scores 0.005552 0.00591 334620188 0.940 0.3472
##
## trial_type = Trial type 4:
## contrast estimate SE df t.ratio p.value
## *D* scores - PI scores 0.003001 0.00611 293120849 0.492 0.6230
##
## Degrees-of-freedom method: kenward-roger
emmeans(fit_disciminability_scoring, list(pairwise ~ DV_type), adjust = "bonferroni")## $`emmeans of DV_type`
## DV_type emmean SE df lower.CL upper.CL
## *D* scores 0.254 0.0109 30357 0.230 0.279
## PI scores 0.250 0.0109 30347 0.226 0.275
##
## Results are averaged over the levels of: trial_type
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
##
## $`pairwise differences of DV_type`
## contrast estimate SE df t.ratio p.value
## *D* scores - PI scores 0.00385 0.00298 320272315 1.289 0.1973
##
## Results are averaged over the levels of: trial_type
## Degrees-of-freedom method: kenward-roger
emmeans(fit_disciminability_scoring, list(pairwise ~ trial_type | DV_type), adjust = "bonferroni")## $`emmeans of trial_type | DV_type`
## DV_type = *D* scores:
## trial_type emmean SE df lower.CL upper.CL
## Trial type 1 0.257 0.0147 84755 0.220 0.294
## Trial type 2 0.236 0.0132 72585 0.203 0.269
## Trial type 3 0.267 0.0174 41481 0.223 0.310
## Trial type 4 0.257 0.0129 104739 0.224 0.289
##
## DV_type = PI scores:
## trial_type emmean SE df lower.CL upper.CL
## Trial type 1 0.257 0.0148 85049 0.220 0.294
## Trial type 2 0.230 0.0132 72895 0.197 0.263
## Trial type 3 0.261 0.0175 41607 0.217 0.305
## Trial type 4 0.254 0.0129 104745 0.221 0.286
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
##
## $`pairwise differences of trial_type | DV_type`
## DV_type = *D* scores:
## contrast estimate SE df t.ratio p.value
## Trial type 1 - Trial type 2 0.020896 0.0159 459177 1.312 1.0000
## Trial type 1 - Trial type 3 -0.009305 0.0180 219859 -0.516 1.0000
## Trial type 1 - Trial type 4 0.000598 0.0159 504342 0.038 1.0000
## Trial type 2 - Trial type 3 -0.030201 0.0145 347025 -2.081 0.2246
## Trial type 2 - Trial type 4 -0.020298 0.0148 682089 -1.375 1.0000
## Trial type 3 - Trial type 4 0.009903 0.0164 218378 0.602 1.0000
##
## DV_type = PI scores:
## contrast estimate SE df t.ratio p.value
## Trial type 1 - Trial type 2 0.027430 0.0160 458847 1.715 0.5184
## Trial type 1 - Trial type 3 -0.003901 0.0181 220162 -0.215 1.0000
## Trial type 1 - Trial type 4 0.003450 0.0159 502959 0.217 1.0000
## Trial type 2 - Trial type 3 -0.031331 0.0146 347919 -2.146 0.1913
## Trial type 2 - Trial type 4 -0.023980 0.0148 681744 -1.621 0.6298
## Trial type 3 - Trial type 4 0.007352 0.0165 218012 0.446 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: bonferroni method for 6 tests
# plot
p_prop_discriminable_scoring <-
ggplot(results_disciminability_scoring, aes(fct_rev(trial_type), estimate, color = DV_type, shape = DV_type, group = DV_type)) +
geom_linerange(aes(ymin = cr_lower, ymax = cr_upper), position = position_dodge(width = dodge_width), linetype = "dotted") +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1, position = position_dodge(width = dodge_width)) +
# geom_line(position = position_dodge(width = dodge_width)) +
geom_point(position = position_dodge(width = dodge_width), size = 2.5) +
scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Worse)", "0.25", "0.50", "0.75", "1.00<br/>(Better)")) +
scale_shape_discrete() +
scale_color_viridis_d(begin = 0.3, end = 0.7) +
mdthemes::md_theme_linedraw() +
scale_x_discrete(labels = c("tt1" = "Trial type 1", "tt2" = "Trial type 2", "tt3" = "Trial type 3", "tt4" = "Trial type 4")) +
labs(x = "",
y = "Proportion of participants who whose scores<br/>differ from one another in pairwise comparisons<br/>",
color = "Scoring method",
shape = "Scoring method") +
theme(legend.position = "top") +
coord_flip(ylim = c(0, 1))
p_prop_discriminable_scoring# fit meta analytic model
fit_disciminability_method <-
lmer(proportion_discriminable ~ 1 + trial_type * method + (method | domain), # non convergence with trial_type as random slope, so had to remove
weights = 1/variance,
data = data_discriminability_D)
# results table
tab_model(fit_disciminability_method,
string.p = "p (corrected)",
ci.hyphen = ", ",
#emph.p = FALSE,
p.adjust = "holm")| proportion discriminable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p (corrected) |
| (Intercept) | 0.26 | 0.23, 0.28 | <0.001 |
| trial_type [tt2] | -0.03 | -0.05, -0.00 | 0.286 |
| trial_type [tt3] | 0.00 | -0.02, 0.03 | 1.000 |
| trial_type [tt4] | -0.01 | -0.03, 0.01 | 1.000 |
| method [basic] | -0.00 | -0.02, 0.02 | 1.000 |
| method [normal] | -0.00 | -0.02, 0.02 | 1.000 |
| method [percent] | -0.00 | -0.02, 0.02 | 1.000 |
|
trial_type [tt2] * method [basic] |
-0.00 | -0.03, 0.03 | 1.000 |
|
trial_type [tt3] * method [basic] |
-0.00 | -0.03, 0.03 | 1.000 |
|
trial_type [tt4] * method [basic] |
-0.00 | -0.03, 0.03 | 1.000 |
|
trial_type [tt2] * method [normal] |
-0.00 | -0.03, 0.03 | 1.000 |
|
trial_type [tt3] * method [normal] |
-0.00 | -0.03, 0.03 | 1.000 |
|
trial_type [tt4] * method [normal] |
-0.00 | -0.03, 0.03 | 1.000 |
|
trial_type [tt2] * method [percent] |
-0.00 | -0.03, 0.03 | 1.000 |
|
trial_type [tt3] * method [percent] |
-0.00 | -0.03, 0.03 | 1.000 |
|
trial_type [tt4] * method [percent] |
-0.00 | -0.03, 0.03 | 1.000 |
| Random Effects | |||
| σ2 | 3.70 | ||
| τ00 domain | 0.00 | ||
| τ11 domain.methodbasic | 0.00 | ||
| τ11 domain.methodnormal | 0.00 | ||
| τ11 domain.methodpercent | 0.00 | ||
| ρ01 | 1.00 | ||
| 0.98 | |||
| 1.00 | |||
| N domain | 35 | ||
| Observations | 560 | ||
| Marginal R2 / Conditional R2 | 0.000 / NA | ||
# summary
summary(fit_disciminability_method)## Linear mixed model fit by REML ['lmerMod']
## Formula: proportion_discriminable ~ 1 + trial_type * method + (method |
## domain)
## Data: data_discriminability_D
## Weights: 1/variance
##
## REML criterion at convergence: -1463.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.83874 -0.57509 0.05272 0.71104 2.34296
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## domain (Intercept) 0.00493695932 0.0702635
## methodbasic 0.00000013075 0.0003616 1.00
## methodnormal 0.00000005142 0.0002268 0.98 0.99
## methodpercent 0.00000011709 0.0003422 1.00 1.00 0.99
## Residual 3.69796864308 1.9230103
## Number of obs: 560, groups: domain, 35
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.25527306 0.01453495 17.563
## trial_typett2 -0.02699398 0.01151785 -2.344
## trial_typett3 0.00388519 0.01160085 0.335
## trial_typett4 -0.00885239 0.01189407 -0.744
## methodbasic -0.00039522 0.01171982 -0.034
## methodnormal -0.00002344 0.01169346 -0.002
## methodpercent -0.00039411 0.01171981 -0.034
## trial_typett2:methodbasic -0.00009045 0.01625505 -0.006
## trial_typett3:methodbasic -0.00056295 0.01637214 -0.034
## trial_typett4:methodbasic -0.00012008 0.01676964 -0.007
## trial_typett2:methodnormal -0.00054747 0.01623850 -0.034
## trial_typett3:methodnormal -0.00123340 0.01636735 -0.075
## trial_typett4:methodnormal -0.00051009 0.01673644 -0.030
## trial_typett2:methodpercent -0.00009243 0.01625505 -0.006
## trial_typett3:methodpercent -0.00056419 0.01637213 -0.034
## trial_typett4:methodpercent -0.00012157 0.01676964 -0.007
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
# plot re
plot_model(fit_disciminability_method,
colors = "bw",
type = "re") +
mdthemes::md_theme_linedraw() +
ggtitle("")# extract re Tau
results_re_tau_disciminability_method <- fit_disciminability_method %>%
merTools::REsdExtract() %>%
#as_tibble(rownames = "trial_type") %>%
as_tibble(rownames = "method") %>%
rename(tau = value) %>%
# mutate(trial_type = case_when(trial_type == "domain.(Intercept)" ~ "Trial type 1",
# trial_type == "domain.trial_typeTrial type 2" ~ "Trial type 2",
# trial_type == "domain.trial_typeTrial type 3" ~ "Trial type 3",
# trial_type == "domain.trial_typeTrial type 4" ~ "Trial type 4"))
mutate(method = case_when(method == "domain.(Intercept)" ~ "bca",
method == "domain.methodbasic" ~ "basic",
method == "domain.methodnormal" ~ "normal",
method == "domain.methodpercent" ~ "percent"))
# extract marginal means
results_emm_disciminability_method <-
summary(emmeans(fit_disciminability_method, ~ method | trial_type)) %>%
dplyr::select(method, trial_type, estimate = emmean, se = SE, ci_lower = lower.CL, ci_upper = upper.CL)
#mutate(trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4"))
# combine
results_disciminability_method <- results_emm_disciminability_method %>%
#full_join(results_re_tau_disciminability_method, by = "trial_type") %>%
full_join(results_re_tau_disciminability_method, by = "method") %>%
mutate(cr_lower = estimate - (1.96 * sqrt(se^2 + tau^2)), # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
cr_upper = estimate + (1.96 * sqrt(se^2 + tau^2))) %>%
mutate(method = fct_relevel(method, "bca", "basic", "normal", "percent"))
# save to disk
write_rds(results_disciminability_method, "models/results_disciminability_method.rds")
# tests
emmeans(fit_disciminability_method, list(pairwise ~ method | trial_type), adjust = "bonferroni")## $`emmeans of method | trial_type`
## trial_type = tt1:
## method emmean SE df lower.CL upper.CL
## bca 0.255 0.0146 149861 0.219 0.292
## basic 0.255 0.0146 149253 0.218 0.291
## normal 0.255 0.0146 148896 0.219 0.292
## percent 0.255 0.0146 149307 0.218 0.291
##
## trial_type = tt2:
## method emmean SE df lower.CL upper.CL
## bca 0.228 0.0144 144363 0.192 0.264
## basic 0.228 0.0144 142730 0.192 0.264
## normal 0.228 0.0144 143139 0.192 0.264
## percent 0.228 0.0144 142778 0.192 0.264
##
## trial_type = tt3:
## method emmean SE df lower.CL upper.CL
## bca 0.259 0.0145 146641 0.223 0.295
## basic 0.258 0.0145 145265 0.222 0.294
## normal 0.258 0.0145 146111 0.222 0.294
## percent 0.258 0.0145 145316 0.222 0.294
##
## trial_type = tt4:
## method emmean SE df lower.CL upper.CL
## bca 0.246 0.0147 154120 0.210 0.283
## basic 0.246 0.0147 151998 0.209 0.283
## normal 0.246 0.0147 151786 0.209 0.283
## percent 0.246 0.0147 152054 0.209 0.283
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
##
## $`pairwise differences of method | trial_type`
## trial_type = tt1:
## contrast estimate SE df t.ratio p.value
## bca - basic 0.000395222 0.0117 754088397 0.034 1.0000
## bca - normal 0.000023437 0.0117 763003410 0.002 1.0000
## bca - percent 0.000394112 0.0117 754435997 0.034 1.0000
## basic - normal -0.000371785 0.0117 760211764 -0.032 1.0000
## basic - percent -0.000001111 0.0117 754028552 0.000 1.0000
## normal - percent 0.000370675 0.0117 760350422 0.032 1.0000
##
## trial_type = tt2:
## contrast estimate SE df t.ratio p.value
## bca - basic 0.000485673 0.0113 883780393 0.043 1.0000
## bca - normal 0.000570911 0.0113 885033942 0.051 1.0000
## bca - percent 0.000486546 0.0113 884214417 0.043 1.0000
## basic - normal 0.000085238 0.0112 896183865 0.008 1.0000
## basic - percent 0.000000873 0.0112 897925167 0.000 1.0000
## normal - percent -0.000084365 0.0112 896353671 -0.008 1.0000
##
## trial_type = tt3:
## contrast estimate SE df t.ratio p.value
## bca - basic 0.000958171 0.0114 832909133 0.084 1.0000
## bca - normal 0.001256840 0.0115 829123868 0.110 1.0000
## bca - percent 0.000958299 0.0114 833299295 0.084 1.0000
## basic - normal 0.000298670 0.0114 836214458 0.026 1.0000
## basic - percent 0.000000128 0.0114 842635860 0.000 1.0000
## normal - percent -0.000298542 0.0114 836348733 -0.026 1.0000
##
## trial_type = tt4:
## contrast estimate SE df t.ratio p.value
## bca - basic 0.000515307 0.0120 688164428 0.043 1.0000
## bca - normal 0.000533523 0.0120 694332865 0.045 1.0000
## bca - percent 0.000515679 0.0120 688416784 0.043 1.0000
## basic - normal 0.000018216 0.0119 704845742 0.002 1.0000
## basic - percent 0.000000372 0.0120 700406298 0.000 1.0000
## normal - percent -0.000017844 0.0119 704958537 -0.001 1.0000
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: bonferroni method for 6 tests
# plot
p_prop_discriminable_method <-
ggplot(results_disciminability_method, aes(trial_type, estimate, color = method, shape = method, group = method)) +
geom_linerange(aes(ymin = cr_lower, ymax = cr_upper), position = position_dodge(width = dodge_width), linetype = "dotted") +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1, position = position_dodge(width = dodge_width)) +
# geom_line(position = position_dodge(width = dodge_width)) +
geom_point(position = position_dodge(width = dodge_width), size = 2.5) +
scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Worse)", "0.25", "0.50", "0.75", "1.00<br/>(Better)")) +
scale_shape_discrete(labels = c("BCA", "Basic", "Normal", "Percentile")) +
scale_color_viridis_d(alpha = 1, begin = 0.3, end = 0.7, labels = c("BCA", "Basic", "Normal", "Percentile")) +
mdthemes::md_theme_linedraw() +
scale_x_discrete(labels = c("tt1" = "Trial type 1", "tt2" = "Trial type 2", "tt3" = "Trial type 3", "tt4" = "Trial type 4")) +
labs(x = "",
y = "Proportion of participants who whose scores<br/>differ from one another in pairwise comparisons",
color = "Bootstrapping method",
shape = "Bootstrapping method") +
theme(legend.position = "top") +
coord_flip(ylim = c(0, 1))
p_prop_discriminable_method## calculate observed ranges
observed_range_estimates_D <- data_estimates_D %>%
group_by(domain, method, trial_type) %>%
dplyr::summarize(min = min(estimate, na.rm = TRUE),
max = max(estimate, na.rm = TRUE),
.groups = "drop") %>%
mutate(range = max - min)
observed_range_estimates_PI <- data_estimates_PI %>%
group_by(domain, method, trial_type) %>%
dplyr::summarize(min = min(estimate, na.rm = TRUE),
max = max(estimate, na.rm = TRUE),
.groups = "drop") %>%
mutate(range = max - min)
# calculate CI / range
data_ci_width_proportions_D <- data_estimates_D %>%
# join this data into the original data
full_join(observed_range_estimates_D, by = c("domain", "method", "trial_type")) %>%
# calculate ci width as a proportion of observed range
mutate(ci_width_proportion = ci_width / range) %>%
mutate(DV_type = "*D* scores")
data_ci_width_proportions_PI <- data_estimates_PI %>%
# join this data into the original data
full_join(observed_range_estimates_PI, by = c("domain", "method", "trial_type")) %>%
# calculate ci width as a proportion of observed range
mutate(ci_width_proportion = ci_width / range) %>%
mutate(DV_type = "PI scores")
# combine
data_ci_width_proportions_combined <-
bind_rows(data_ci_width_proportions_D,
data_ci_width_proportions_PI) %>%
mutate(domain = as.factor(domain),
method = fct_relevel(method, "bca", "basic", "normal", "percent"),
trial_type = fct_relevel(trial_type, "tt1", "tt2", "tt3", "tt4"))# fit model
fit_ci_width_proportions_scoring <-
data_ci_width_proportions_combined %>%
filter(method == "bca" & DV_type %in% c("*D* scores", "PI scores")) %>%
mutate(trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
trial_type == "tt2" ~ "Trial type 2",
trial_type == "tt3" ~ "Trial type 3",
trial_type == "tt4" ~ "Trial type 4")) %>%
lmer(ci_width_proportion ~ trial_type * DV_type + (trial_type | domain) + (1 | unique_id),
data = .)
# results table
tab_model(fit_ci_width_proportions_scoring,
string.p = "p (corrected)",
ci.hyphen = ", ",
#emph.p = FALSE,
p.adjust = "holm")| ci width proportion | |||
|---|---|---|---|
| Predictors | Estimates | CI | p (corrected) |
| (Intercept) | 0.83 | 0.76, 0.90 | <0.001 |
| trial_type [Trial type 2] | 0.05 | -0.01, 0.12 | 0.416 |
| trial_type [Trial type 3] | 0.00 | -0.09, 0.09 | 1.000 |
| trial_type [Trial type 4] | -0.03 | -0.10, 0.05 | 1.000 |
| DV_type PI scores | -0.02 | -0.03, -0.02 | <0.001 |
|
trial_type [Trial type 2] * DV_type PI scores |
-0.00 | -0.01, 0.00 | 1.000 |
|
trial_type [Trial type 3] * DV_type PI scores |
-0.02 | -0.03, -0.01 | 0.002 |
|
trial_type [Trial type 4] * DV_type PI scores |
-0.00 | -0.01, 0.01 | 1.000 |
| Random Effects | |||
| σ2 | 0.01 | ||
| τ00 unique_id | 0.00 | ||
| τ00 domain | 0.05 | ||
| τ11 domain.trial_typeTrial type 2 | 0.03 | ||
| τ11 domain.trial_typeTrial type 3 | 0.07 | ||
| τ11 domain.trial_typeTrial type 4 | 0.05 | ||
| ρ01 domain.trial_typeTrial type 2 | -0.49 | ||
| ρ01 domain.trial_typeTrial type 3 | -0.47 | ||
| ρ01 domain.trial_typeTrial type 4 | -0.71 | ||
| ICC | 0.87 | ||
| N domain | 35 | ||
| N unique_id | 1464 | ||
| Observations | 11712 | ||
| Marginal R2 / Conditional R2 | 0.018 / 0.870 | ||
# # plot RE effects
# plot_model(fit_ci_width_proportions_scoring,
# colors = "bw",
# type = "re")[[1]] +
# mdthemes::md_theme_linedraw() +
# ggtitle("")
# extract re Tau
results_re_tau_ci_width_proportions_scoring <-
merTools::REsdExtract(fit_ci_width_proportions_scoring) %>%
as_tibble(rownames = "trial_type") %>%
rename(tau = value) %>%
mutate(trial_type = case_when(trial_type == "domain.(Intercept)" ~ "Trial type 1",
trial_type == "domain.trial_typeTrial type 2" ~ "Trial type 2",
trial_type == "domain.trial_typeTrial type 3" ~ "Trial type 3",
trial_type == "domain.trial_typeTrial type 4" ~ "Trial type 4")) %>%
drop_na()
# extract marginal means
results_emm_ci_width_proportions_scoring <-
summary(emmeans(fit_ci_width_proportions_scoring, ~ DV_type | trial_type)) %>%
dplyr::select(DV_type, trial_type, estimate = emmean, se = SE, ci_lower = asymp.LCL, ci_upper = asymp.UCL) %>%
mutate(DV_type = fct_relevel(DV_type, "*D* scores", "PI scores"),
trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4"))
# combine
results_ci_width_proportions_scoring <- results_emm_ci_width_proportions_scoring %>%
full_join(results_re_tau_ci_width_proportions_scoring, by = "trial_type") %>%
mutate(cr_lower = estimate - (1.96 * sqrt(se^2 + tau^2)), # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
cr_upper = estimate + (1.96 * sqrt(se^2 + tau^2)))
# save to disk
write_rds(results_ci_width_proportions_scoring, "models/results_ci_width_proportions_scoring.rds")
# tests
emmeans(fit_ci_width_proportions_scoring, list(pairwise ~ DV_type | trial_type), adjust = "bonferroni")## $`emmeans of DV_type | trial_type`
## trial_type = Trial type 1:
## DV_type emmean SE df asymp.LCL asymp.UCL
## *D* scores 0.830 0.0382 Inf 0.745 0.916
## PI scores 0.808 0.0382 Inf 0.723 0.894
##
## trial_type = Trial type 2:
## DV_type emmean SE df asymp.LCL asymp.UCL
## *D* scores 0.885 0.0357 Inf 0.805 0.965
## PI scores 0.858 0.0357 Inf 0.778 0.938
##
## trial_type = Trial type 3:
## DV_type emmean SE df asymp.LCL asymp.UCL
## *D* scores 0.831 0.0440 Inf 0.732 0.930
## PI scores 0.793 0.0440 Inf 0.694 0.891
##
## trial_type = Trial type 4:
## DV_type emmean SE df asymp.LCL asymp.UCL
## *D* scores 0.804 0.0287 Inf 0.740 0.869
## PI scores 0.782 0.0287 Inf 0.718 0.846
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 2 estimates
##
## $`pairwise differences of DV_type | trial_type`
## trial_type = Trial type 1:
## contrast estimate SE df z.ratio p.value
## *D* scores - PI scores 0.0220 0.00323 Inf 6.805 <.0001
##
## trial_type = Trial type 2:
## contrast estimate SE df z.ratio p.value
## *D* scores - PI scores 0.0267 0.00323 Inf 8.252 <.0001
##
## trial_type = Trial type 3:
## contrast estimate SE df z.ratio p.value
## *D* scores - PI scores 0.0384 0.00323 Inf 11.883 <.0001
##
## trial_type = Trial type 4:
## contrast estimate SE df z.ratio p.value
## *D* scores - PI scores 0.0221 0.00323 Inf 6.852 <.0001
##
## Degrees-of-freedom method: asymptotic
emmeans(fit_ci_width_proportions_scoring, list(pairwise ~ trial_type | DV_type), adjust = "bonferroni")## $`emmeans of trial_type | DV_type`
## DV_type = *D* scores:
## trial_type emmean SE df asymp.LCL asymp.UCL
## Trial type 1 0.830 0.0382 Inf 0.735 0.925
## Trial type 2 0.885 0.0357 Inf 0.796 0.974
## Trial type 3 0.831 0.0440 Inf 0.721 0.941
## Trial type 4 0.804 0.0287 Inf 0.733 0.876
##
## DV_type = PI scores:
## trial_type emmean SE df asymp.LCL asymp.UCL
## Trial type 1 0.808 0.0382 Inf 0.713 0.903
## Trial type 2 0.858 0.0357 Inf 0.769 0.947
## Trial type 3 0.793 0.0440 Inf 0.683 0.902
## Trial type 4 0.782 0.0287 Inf 0.710 0.854
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
##
## $`pairwise differences of trial_type | DV_type`
## DV_type = *D* scores:
## contrast estimate SE df z.ratio p.value
## Trial type 1 - Trial type 2 -0.054947 0.0317 Inf -1.733 0.4991
## Trial type 1 - Trial type 3 -0.000945 0.0460 Inf -0.021 1.0000
## Trial type 1 - Trial type 4 0.025783 0.0374 Inf 0.690 1.0000
## Trial type 2 - Trial type 3 0.054002 0.0370 Inf 1.459 0.8676
## Trial type 2 - Trial type 4 0.080730 0.0306 Inf 2.642 0.0495
## Trial type 3 - Trial type 4 0.026728 0.0366 Inf 0.730 1.0000
##
## DV_type = PI scores:
## contrast estimate SE df z.ratio p.value
## Trial type 1 - Trial type 2 -0.050271 0.0317 Inf -1.585 0.6777
## Trial type 1 - Trial type 3 0.015466 0.0460 Inf 0.336 1.0000
## Trial type 1 - Trial type 4 0.025937 0.0374 Inf 0.694 1.0000
## Trial type 2 - Trial type 3 0.065738 0.0370 Inf 1.776 0.4545
## Trial type 2 - Trial type 4 0.076208 0.0306 Inf 2.494 0.0759
## Trial type 3 - Trial type 4 0.010470 0.0366 Inf 0.286 1.0000
##
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 6 tests
# plot
dodge_width <- 0.8
p_ci_width_proportion_observed_range_scoring <-
ggplot(results_ci_width_proportions_scoring, aes(fct_rev(trial_type), estimate, color = DV_type, shape = DV_type, group = DV_type)) +
geom_linerange(aes(ymin = cr_lower, ymax = cr_upper), position = position_dodge(width = dodge_width), linetype = "dotted") +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1, position = position_dodge(width = dodge_width)) +
# geom_line(position = position_dodge(width = dodge_width)) +
geom_point(position = position_dodge(width = dodge_width), size = 2.5) +
scale_shape_discrete(labels = c("*D* scores", "PI scores")) +
scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Better)", "0.25", "0.50", "0.75", "1.00<br/>(Worse)")) +
scale_color_viridis_d(begin = 0.3, end = 0.7, labels = c("*D* scores", "PI scores")) +
mdthemes::md_theme_linedraw() +
scale_x_discrete(labels = c("tt1" = "Trial type 1", "tt2" = "Trial type 2", "tt3" = "Trial type 3", "tt4" = "Trial type 4")) +
labs(x = "",
y = "Within-subject 95% CI widths /<br/>between-subjects observed range<br/>",
color = "Scoring method",
shape = "Scoring method") +
theme(legend.position = "top") +
coord_flip(ylim = c(0, 1))
p_ci_width_proportion_observed_range_scoringdata_ci_width_proportions_combined_method <-
data_ci_width_proportions_D %>%
mutate(domain = as.factor(domain),
method = as.factor(method),
method = fct_relevel(method, "bca", "basic", "normal", "percent"),
trial_type = fct_relevel(trial_type, "tt1", "tt2", "tt3", "tt4")) %>%
select(unique_id, domain, trial_type, method, ci_width_proportion)
# fit model
if(file.exists("models/fit_ci_width_proportions_method.rds")) {
fit_ci_width_proportions_method <- read_rds("models/fit_ci_width_proportions_method.rds")
} else {
fit_ci_width_proportions_method <-
lmer(ci_width_proportion ~ trial_type * method + (trial_type | domain) + (1 | unique_id),
data = data_ci_width_proportions_combined_method)
write_rds(fit_ci_width_proportions_method, "models/fit_ci_width_proportions_method.rds")
}
# results table
tab_model(fit_ci_width_proportions_method,
string.p = "p (corrected)",
ci.hyphen = ", ",
#emph.p = FALSE,
p.adjust = "holm")| ci width proportion | |||
|---|---|---|---|
| Predictors | Estimates | CI | p (corrected) |
| (Intercept) | 0.83 | 0.75, 0.90 | <0.001 |
| trial_type [tt2] | 0.06 | -0.01, 0.12 | 1.000 |
| trial_type [tt3] | 0.00 | -0.10, 0.11 | 1.000 |
| trial_type [tt4] | -0.02 | -0.10, 0.06 | 1.000 |
| method [basic] | -0.01 | -0.02, -0.01 | 0.001 |
| method [normal] | -0.01 | -0.02, -0.00 | 0.051 |
| method [percent] | -0.01 | -0.02, -0.01 | 0.001 |
|
trial_type [tt2] * method [basic] |
0.01 | -0.00, 0.02 | 1.000 |
|
trial_type [tt3] * method [basic] |
0.01 | -0.00, 0.01 | 1.000 |
|
trial_type [tt4] * method [basic] |
0.01 | -0.00, 0.02 | 1.000 |
|
trial_type [tt2] * method [normal] |
0.01 | -0.00, 0.01 | 1.000 |
|
trial_type [tt3] * method [normal] |
0.00 | -0.00, 0.01 | 1.000 |
|
trial_type [tt4] * method [normal] |
0.01 | -0.00, 0.02 | 1.000 |
|
trial_type [tt2] * method [percent] |
0.01 | -0.00, 0.02 | 1.000 |
|
trial_type [tt3] * method [percent] |
0.01 | -0.00, 0.01 | 1.000 |
|
trial_type [tt4] * method [percent] |
0.01 | -0.00, 0.02 | 1.000 |
| Random Effects | |||
| σ2 | 0.01 | ||
| τ00 unique_id | 0.00 | ||
| τ00 domain | 0.05 | ||
| τ11 domain.trial_typett2 | 0.04 | ||
| τ11 domain.trial_typett3 | 0.09 | ||
| τ11 domain.trial_typett4 | 0.06 | ||
| ρ01 domain.trial_typett2 | -0.52 | ||
| ρ01 domain.trial_typett3 | -0.54 | ||
| ρ01 domain.trial_typett4 | -0.68 | ||
| ICC | 0.87 | ||
| N domain | 35 | ||
| N unique_id | 1464 | ||
| Observations | 23424 | ||
| Marginal R2 / Conditional R2 | 0.013 / 0.876 | ||
# summary
summary(fit_ci_width_proportions_method)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: ci_width_proportion ~ trial_type * method + (trial_type | domain) +
## (1 | unique_id)
## Data: data_ci_width_proportions_combined_method
##
## REML criterion at convergence: -42482.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7931 -0.4792 0.0852 0.5727 5.7306
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## unique_id (Intercept) 0.004573 0.06762
## domain (Intercept) 0.051563 0.22707
## trial_typett2 0.042396 0.20590 -0.52
## trial_typett3 0.093386 0.30559 -0.54 0.59
## trial_typett4 0.058856 0.24260 -0.68 0.53 0.72
## Residual 0.007930 0.08905
## Number of obs: 23424, groups: unique_id, 1464; domain, 35
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 0.827669 0.038526 33.977736 21.483
## trial_typett2 0.056406 0.034981 34.341297 1.612
## trial_typett3 0.004559 0.051774 34.100392 0.088
## trial_typett4 -0.020482 0.041158 34.227455 -0.498
## methodbasic -0.013544 0.003291 21842.663130 -4.115
## methodnormal -0.009487 0.003291 21842.663134 -2.882
## methodpercent -0.013544 0.003291 21842.663133 -4.115
## trial_typett2:methodbasic 0.006058 0.004655 21842.663131 1.301
## trial_typett3:methodbasic 0.005098 0.004655 21842.663123 1.095
## trial_typett4:methodbasic 0.006675 0.004655 21842.663124 1.434
## trial_typett2:methodnormal 0.005841 0.004655 21842.663133 1.255
## trial_typett3:methodnormal 0.004693 0.004655 21842.663125 1.008
## trial_typett4:methodnormal 0.006251 0.004655 21842.663128 1.343
## trial_typett2:methodpercent 0.006058 0.004655 21842.663132 1.301
## trial_typett3:methodpercent 0.005098 0.004655 21842.663126 1.095
## trial_typett4:methodpercent 0.006675 0.004655 21842.663127 1.434
## Pr(>|t|)
## (Intercept) < 0.0000000000000002 ***
## trial_typett2 0.11601
## trial_typett3 0.93035
## trial_typett4 0.62192
## methodbasic 0.0000389 ***
## methodnormal 0.00395 **
## methodpercent 0.0000389 ***
## trial_typett2:methodbasic 0.19312
## trial_typett3:methodbasic 0.27345
## trial_typett4:methodbasic 0.15158
## trial_typett2:methodnormal 0.20954
## trial_typett3:methodnormal 0.31336
## trial_typett4:methodnormal 0.17934
## trial_typett2:methodpercent 0.19312
## trial_typett3:methodpercent 0.27345
## trial_typett4:methodpercent 0.15158
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00321383 (tol = 0.002, component 1)
# # plot re
# plot_model(fit_ci_width_proportions_method,
# colors = "bw",
# type = "re") +
# mdthemes::md_theme_linedraw() +
# ggtitle("")
# extract re Tau
results_re_tau_ci_width_proportions_method <- fit_ci_width_proportions_method %>%
merTools::REsdExtract() %>%
as_tibble(rownames = "trial_type") %>%
rename(tau = value) %>%
mutate(trial_type = case_when(trial_type == "domain.(Intercept)" ~ "Trial type 1",
trial_type == "domain.trial_typett2" ~ "Trial type 2",
trial_type == "domain.trial_typett3" ~ "Trial type 3",
trial_type == "domain.trial_typett4" ~ "Trial type 4"))
# extract marginal means
results_emm_ci_width_proportions_method <-
summary(emmeans(fit_ci_width_proportions_method, ~ method | trial_type)) %>%
dplyr::select(method, trial_type, estimate = emmean, se = SE, ci_lower = asymp.LCL, ci_upper = asymp.UCL) %>%
mutate(trial_type = case_when(trial_type == "tt1" ~ "Trial type 1",
trial_type == "tt2" ~ "Trial type 2",
trial_type == "tt3" ~ "Trial type 3",
trial_type == "tt4" ~ "Trial type 4"),
trial_type = fct_relevel(trial_type, "Trial type 1", "Trial type 2", "Trial type 3", "Trial type 4"))
# combine
results_ci_width_proportions_method <- results_emm_ci_width_proportions_method %>%
full_join(results_re_tau_ci_width_proportions_method, by = "trial_type") %>%
mutate(cr_lower = estimate - (1.96 * sqrt(se^2 + tau^2)), # as in metafor package's implementation of credibility intervals, see metafor::predict.rma.R
cr_upper = estimate + (1.96 * sqrt(se^2 + tau^2))) %>%
drop_na()
# save to disk
write_rds(results_ci_width_proportions_method, "models/results_ci_width_proportions_method.rds")
# tests
emmeans(fit_ci_width_proportions_method, list(pairwise ~ method | trial_type), adjust = "bonferroni")## $`emmeans of method | trial_type`
## trial_type = tt1:
## method emmean SE df asymp.LCL asymp.UCL
## bca 0.828 0.0385 Inf 0.731 0.924
## basic 0.814 0.0385 Inf 0.718 0.910
## normal 0.818 0.0385 Inf 0.722 0.914
## percent 0.814 0.0385 Inf 0.718 0.910
##
## trial_type = tt2:
## method emmean SE df asymp.LCL asymp.UCL
## bca 0.884 0.0362 Inf 0.794 0.974
## basic 0.877 0.0362 Inf 0.786 0.967
## normal 0.880 0.0362 Inf 0.790 0.971
## percent 0.877 0.0362 Inf 0.786 0.967
##
## trial_type = tt3:
## method emmean SE df asymp.LCL asymp.UCL
## bca 0.832 0.0449 Inf 0.720 0.944
## basic 0.824 0.0449 Inf 0.712 0.936
## normal 0.827 0.0449 Inf 0.715 0.940
## percent 0.824 0.0449 Inf 0.712 0.936
##
## trial_type = tt4:
## method emmean SE df asymp.LCL asymp.UCL
## bca 0.807 0.0320 Inf 0.727 0.887
## basic 0.800 0.0320 Inf 0.720 0.880
## normal 0.804 0.0320 Inf 0.724 0.884
## percent 0.800 0.0320 Inf 0.720 0.880
##
## Degrees-of-freedom method: asymptotic
## Confidence level used: 0.95
## Conf-level adjustment: bonferroni method for 4 estimates
##
## $`pairwise differences of method | trial_type`
## trial_type = tt1:
## contrast estimate SE df z.ratio p.value
## bca - basic 0.01354 0.00329 Inf 4.115 0.0002
## bca - normal 0.00949 0.00329 Inf 2.882 0.0237
## bca - percent 0.01354 0.00329 Inf 4.115 0.0002
## basic - normal -0.00406 0.00329 Inf -1.233 1.0000
## basic - percent 0.00000 0.00329 Inf 0.000 1.0000
## normal - percent 0.00406 0.00329 Inf 1.233 1.0000
##
## trial_type = tt2:
## contrast estimate SE df z.ratio p.value
## bca - basic 0.00749 0.00329 Inf 2.274 0.1376
## bca - normal 0.00365 0.00329 Inf 1.108 1.0000
## bca - percent 0.00749 0.00329 Inf 2.274 0.1376
## basic - normal -0.00384 0.00329 Inf -1.167 1.0000
## basic - percent 0.00000 0.00329 Inf 0.000 1.0000
## normal - percent 0.00384 0.00329 Inf 1.167 1.0000
##
## trial_type = tt3:
## contrast estimate SE df z.ratio p.value
## bca - basic 0.00845 0.00329 Inf 2.566 0.0617
## bca - normal 0.00479 0.00329 Inf 1.456 0.8717
## bca - percent 0.00845 0.00329 Inf 2.566 0.0617
## basic - normal -0.00365 0.00329 Inf -1.110 1.0000
## basic - percent 0.00000 0.00329 Inf 0.000 1.0000
## normal - percent 0.00365 0.00329 Inf 1.110 1.0000
##
## trial_type = tt4:
## contrast estimate SE df z.ratio p.value
## bca - basic 0.00687 0.00329 Inf 2.087 0.2214
## bca - normal 0.00324 0.00329 Inf 0.983 1.0000
## bca - percent 0.00687 0.00329 Inf 2.087 0.2214
## basic - normal -0.00363 0.00329 Inf -1.104 1.0000
## basic - percent 0.00000 0.00329 Inf 0.000 1.0000
## normal - percent 0.00363 0.00329 Inf 1.104 1.0000
##
## Degrees-of-freedom method: asymptotic
## P value adjustment: bonferroni method for 6 tests
# plot
p_ci_width_proportion_observed_range_method <-
ggplot(results_ci_width_proportions_method, aes(trial_type, estimate,
color = method,
shape = method,
group = method)) +
geom_linerange(aes(ymin = cr_lower, ymax = cr_upper), position = position_dodge(width = dodge_width), linetype = "dotted") +
geom_linerange(aes(ymin = ci_lower, ymax = ci_upper), size = 1, position = position_dodge(width = dodge_width)) +
# geom_line(position = position_dodge(width = dodge_width)) +
geom_point(position = position_dodge(width = dodge_width), size = 2.5) +
mdthemes::md_theme_linedraw() +
scale_y_continuous(breaks = c(0, .25, .5, .75, 1), labels = c("0.00<br/>(Better)", "0.25", "0.50", "0.75", "1.00<br/>(Worse)")) +
scale_shape_discrete(labels = c("BCA", "Basic", "Normal", "Percentile")) +
scale_color_viridis_d(alpha = 1, begin = 0.3, end = 0.7, labels = c("BCA", "Basic", "Normal", "Percentile")) +
labs(shape = "Bootstrapping method",
color = "Bootstrapping method",
x = "",
y = "95% CI width / observed range") +
theme(legend.position = "top") +
coord_flip(ylim = c(0, 1))
p_ci_width_proportion_observed_range_methodPlot 1 is merely illusatrative. It shows the bootstrapped CIs for all participants, split by domain, but not splitting by trial type. D scores and the BCA method are displayed as they are used for the subsequent analyses.
p_cis_by_domainggsave(filename = "plots/figure_1_cis_by_domain.pdf",
plot = p_cis_by_domain,
device = "pdf",
# path = NULL,
# dpi = 300,
units = "in",
width = 8,
height = 8,
limitsize = TRUE)Most probable CI width for D scores when bootstrapped using four different methods. Very similar results are found across methods. Overall Maximum A-Posteroi width (MAP, i.e., the mode of a continuous variable) was D = 1.31 +/- 0.01 between bootstrapping methods. Some domains and trial types did demonstrate smaller most probable widths.
NB I elected not to meta-analyze these widths as they demonstrate very large skew at the individual level, which violate the assumptions of linear meta-analysis and underestimate the typical width (ie estimated mean widths << MAP observed widths). Rather than meta analyze, I simply report the domain and trial type level MAP values. More informative and valid analyses are presented below - ones which can directly compare the D and PI as an alternative. That could not be accomplished with a direct comparison with D/PI scores’ 95% CIs as they are on different scales and follow different distributions.
Given the minimal differences between bootstrapping methods, I employ the Bias Corrected and Accelerated methods (BCA) for all subsequent primary analyses. Sensitivity analyses that compare between bootstrapping methods are computed but not included in the main manuscript.
p_ci_widthsggsave(filename = "plots/figure_2_ci_widths.pdf",
plot = p_ci_widths,
device = "pdf",
# path = NULL,
# dpi = 300,
units = "in",
width = 8,
height = 6,
limitsize = TRUE)The results of three hierarchical/meta analytic models are presented below, all of which provide information via different methods regarding how informative an individual participants’ D (or PI) score is in terms of being able to state that they demonstrated evidence of a bias/IRAP effect/implicit attitude, whether that individual can be discriminated from other individuals in the same domain using the same trial type, and how much of the range of observed scores an individuals score’s CI spans.
All meta analyses compare D and PI scores to assess whether the results are dependent on the D algorithm which has been argued to be suboptimal. That is, I assess whether this issue can be trivially resolved by scoring the data a different way.
Note that the theoretical max possible range of D scores is -2 to +2, but such extreme scores are practically impossible. As such, in order to understand the CI width in terms of realistic data - and also in order to compare D and PI which are on different scales and distributions - I standardize CI widths by the observed range of scores for each domain and trial type.
p_combined <-
p_prop_nonzero_scoring +
p_prop_discriminable_scoring +
p_ci_width_proportion_observed_range_scoring +
plot_layout(ncol = 1, guides = "collect") & theme(legend.position = "bottom")
p_combinedggsave(filename = "plots/figure_3_metaanalyses.pdf",
plot = p_combined,
device = "pdf",
# path = NULL,
# dpi = 300,
units = "in",
width = 5,
height = 7,
limitsize = TRUE)sessionInfo()## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_IE.UTF-8/en_IE.UTF-8/en_IE.UTF-8/C/en_IE.UTF-8/en_IE.UTF-8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] lmerTest_3.1-3 ggstance_0.3.4 emmeans_1.4.6 sjPlot_2.8.4
## [5] lme4_1.1-26 Matrix_1.2-18 mdthemes_0.1.0 patchwork_1.0.1
## [9] bayestestR_0.7.5 boot_1.3-27 kableExtra_1.3.1 knitr_1.30
## [13] forcats_0.5.0 stringr_1.4.0 dplyr_1.0.3 purrr_0.3.4
## [17] readr_1.3.1 tidyr_1.1.2 tibble_3.1.0 ggplot2_3.3.3
## [21] tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.0-10 minqa_1.2.4 colorspace_2.0-0
## [4] ellipsis_0.3.1 sjlabelled_1.1.7 estimability_1.3
## [7] snakecase_0.11.0 markdown_1.1 parameters_0.8.6
## [10] fs_1.4.1 gridtext_0.1.1 ggtext_0.1.0.9000
## [13] rstudioapi_0.13 glmmTMB_1.0.2.1 farver_2.1.0
## [16] fansi_0.4.2 mvtnorm_1.1-1 lubridate_1.7.9
## [19] xml2_1.3.2 codetools_0.2-16 splines_4.0.2
## [22] sjmisc_2.8.5 jsonlite_1.7.2 nloptr_1.2.2.2
## [25] ggeffects_0.14.3 pbkrtest_0.4-8.6 broom_0.7.2
## [28] dbplyr_1.4.3 broom.mixed_0.2.6 shiny_1.5.0
## [31] effectsize_0.4.0 compiler_4.0.2 httr_1.4.1
## [34] sjstats_0.18.0 backports_1.1.9 fastmap_1.0.1
## [37] assertthat_0.2.1 cli_2.3.1 later_1.0.0
## [40] htmltools_0.5.0 tools_4.0.2 coda_0.19-3
## [43] gtable_0.3.0 glue_1.4.2 reshape2_1.4.4
## [46] merTools_0.5.2 Rcpp_1.0.6 cellranger_1.1.0
## [49] vctrs_0.3.6 nlme_3.1-148 iterators_1.0.12
## [52] insight_0.10.0 xfun_0.19 rvest_0.3.5
## [55] mime_0.9 lifecycle_1.0.0 statmod_1.4.35
## [58] MASS_7.3-53 zoo_1.8-8 scales_1.1.1
## [61] promises_1.1.0 hms_0.5.3 sandwich_2.5-1
## [64] TMB_1.7.18 yaml_2.2.1 stringi_1.5.3
## [67] highr_0.8 foreach_1.5.0 plotrix_3.7-8
## [70] blme_1.0-5 rlang_0.4.10 pkgconfig_2.0.3
## [73] arm_1.11-1 evaluate_0.14 lattice_0.20-41
## [76] labeling_0.4.2 tidyselect_1.1.0 plyr_1.8.6
## [79] magrittr_2.0.1 R6_2.5.0 generics_0.1.0
## [82] multcomp_1.4-13 DBI_1.1.0 pillar_1.5.1
## [85] haven_2.3.1 withr_2.4.1 abind_1.4-5
## [88] survival_3.1-12 performance_0.4.6 janitor_2.0.1
## [91] modelr_0.1.8 crayon_1.4.1 utf8_1.2.1
## [94] rmarkdown_2.5 grid_4.0.2 readxl_1.3.1
## [97] reprex_0.3.0 digest_0.6.27 webshot_0.5.2
## [100] xtable_1.8-4 numDeriv_2016.8-1.1 httpuv_1.5.2
## [103] munsell_0.5.0 viridisLite_0.3.0